gmac-lab-chart

1 📕 READ ME

SUMMARY
This R code is used to estimate the relationship between oxygen consumption (MO₂) and ambient oxygen partial pressure (PO₂) in the Common Galaxias (Galaxias maculatus). It also estimates the critical partial pressure of oxygen for aerobic metabolism (Pcrit), which is commonly understood as the threshold below which the oxygen consumption rate can no longer be sustained. The associated article is “The role of osmorespiratory compromise in hypoxia tolerance of the purportedly oxyconforming teleost Galaxias maculatus.” If you click the Code button in the top right of the HTML document, you can download the .Rmd file.

AIM
The article aims to test whether Galaxias maculatus can maintain oxygen consumption (MO₂) as ambient PO₂ falls and, if so, at what level it reaches the critical partial pressure of oxygen for aerobic metabolism (Pcrit).

AUTHORS
Timothy D. Clark ^[a]
Luis L. Kuchenmüller ^[a]
Elizabeth C. Hoots ^[a]
Maryane Gradito ^[a]
Jake M. Martin ^[a,b]
🚧 In no particular order 🚧

AFFILIATIONS
[a] School of Life and Environmental Sciences, Deakin University, Geelong, VIC, Australia
[b] School of Biological Sciences, Monash University, Clayton, Victoria, Australia

CONTRIBUTOR ROLES
🚧 To be added 🚧 Based on the Contributor Roles Taxonomy (CRediT)

DISCLAIMER
I (Jake Martin) am dyslexic. I have made an effort to review the script for grammatical errors, but some will likely remain. I apologise. Please reach out via the contact details below if anything is unclear.

2 📧 Contact

Jake M. Martin

📧 Email: jake.martin@deakin.edu.au

📧 Alt Email: jake.martin.research@gmail.com

🌐 Web: jakemartin.org

🐙 GitHub: JakeMartinResearch

3 ⚙️ Knit settings

Here we define our Knit settings, to make the output more user friendly, and to cache output for faster knitting.

#kniter seetting
knitr::opts_chunk$set(
message = FALSE,
warning = FALSE, # no warnings
cache = TRUE,# Cacheing to save time when kniting
tidy = TRUE
)

4 📦 Required packages

These are the R packages required for this script. You will need to install a package called pacman to run the p_load function.

# you need to install pacman first install.packages("pacman")
pacman::p_load("ggplot2", 
               "ggthemes", 
               "ggfortify", 
               "gtExtras", 
               "igraph",
               "dagitty",
               "ggdag",
               "ggridges",
               "gghalves",
               "ggExtra",
               "gridExtra",
               "corrplot",
               "RColorBrewer", 
               "gt", 
               "gtsummary",
               "grid",
               "plotly",
               "bayesplot", # data visualisation
               
                "tidyverse", 
               "janitor", 
               "readxl", 
               "broom.mixed", 
               "data.table", 
               "devtools",
               "hms", # data tidy
               
               "marginaleffects", 
               "brms", 
               "rstan", 
               "performance", 
               "emmeans", 
               "tidybayes", 
               "vegan",
               "betareg",
               "lme4", 
               "car", 
               "lmerTest",
               "qqplotr",
               "respirometry",
               "mclust",
               "furrr", # modelling 
              
               
               "datawizard", 
               "SRS", # data manipulation 
               
               "sessioninfo" # reproducibility
                       )

5 🔧 Custom functions

Here are some custom function used within this script.

bayes_incremental_regression_by_id(): A custom function to build Bayesian incremental regressions. It is formatted for running a list of sub group models (ids) in parallel with 4 cores. It uses brm() with a Gaussian error function.

Use: The function takes an grouping factor/id that will be used to filter the data for the regression, if none is proved it will use all data (id_i), the column name of the grouping factor/charter in the data frame, if none is proved it will use all data (id_name), the data frame (data), the predictor of interest (predictor), the response of interest (response), the random seed number for model reproducibility (seed_number), where you would like to save the model outputs (save_models = logical argument), the output directory for the saved rds models (mod_output_wd)

bayes_incremental_regression_by_id(id_i = “your_id”, id_name = “column name for ids”, data = “your data”, predictor = “predictor variable”, response = “response variable”, seed_number = integer, save_models = TRUE/FALSE, mod_output_wd = “directory for model output”)

# Instead we could also use a distributional regression approach, by
# specifically modelling the variance by DO (e.g. sigma ~ DO). Weighting may
# not be required in this case, I don't think higher density of vaules in a
# given space will effect Bayesian estimates like it does in frequentist
# models. See discourse https://discourse.mc-stan.org/t/weights-in-brm/4278

bayes_incremental_regression_by_id <- function(id_i, id_name, data, predictor, response,
    seed_number, save_models, mod_output_wd) {
    # Initiate an empty list to store models
    models <- list()

    # Check if id_name is missing, NULL, or blank, and assign NA if so
    if (missing(id_name) || is.null(id_name) || id_name == "") {
        id_name <- NA
    }

    # Check if id_i is missing, NULL, or blank, and assign NA if so
    if (missing(id_i) || is.null(id_i) || id_i == "") {
        id_name <- NA
    }

    # Filter data for the current ID if id_name is given as a factor or
    # character and id_i is defined
    df_i <- data %>%
        dplyr::filter(if (!is.na(id_i) && (is.factor(data[[id_name]]) || is.character(data[[id_name]]))) {
            !!rlang::sym(id_name) == id_i
        } else {
            TRUE
        })

    # Dynamically create formulas
    formula_lm_0 <- reformulate("1", response)
    formula_lm_1 <- reformulate(predictor, response)
    formula_lm_2 <- reformulate(sprintf("poly(%s, 2)", predictor), response)
    formula_lm_3 <- reformulate(sprintf("poly(%s, 3)", predictor), response)

    # Fit and store models in the list
    models[[paste0(id_i, "_lm_0")]] <- brm(bf(formula_lm_0, family = gaussian()),
        data = df_i, cores = 4, seed = seed_number, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_0"))

    models[[paste0(id_i, "_lm_1")]] <- brm(bf(formula_lm_1, family = gaussian()),
        data = df_i, cores = 4, seed = seed_number, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_1"))

    models[[paste0(id_i, "_lm_2")]] <- brm(bf(formula_lm_2, family = gaussian()),
        data = df_i, cores = 4, seed = seed_number, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_2"))

    models[[paste0(id_i, "_lm_3")]] <- brm(bf(formula_lm_3, family = gaussian()),
        data = df_i, cores = 4, seed = 143019, save_pars = save_pars(all = save_models),
        sample_prior = FALSE, silent = TRUE, file = paste0(mod_output_wd, "/", id_i,
            "_lm_3"))

    # Return the list of models for the current ID
    return(models)
}

load_rds9(): A custom function to load all rds models in a directory and store in a list. Takes a directory with .rds files

load_rds <- function(model_dw) {
    # List all .rds files in the directory
    model_file_list <- list.files(path = model_dw, pattern = "\\.rds$", full.names = TRUE)

    # Initialise an empty list to store models
    model_store_list <- list()

    # Iterate through each file and load the RDS
    for (mod_i in model_file_list) {
        mod <- readRDS(file = mod_i)  # Read the RDS file
        model_name <- tools::file_path_sans_ext(basename(mod_i))  # Extract the file name without extension
        model_store_list[[model_name]] <- mod  # Store it in the list
    }

    # Return the list of models
    return(model_store_list)
}

incremental_regression_bayes_fits(): A custom function for pulling model fits, loo and r2 using loo() and bayes_R2(), respectively. Takes a list of brm models.

# Define Function to Process the data for each ID
incremental_regression_bayes_fits <- function(models) {
  
  loo_results_list <- list()
  
  # Iterate over the names of the models
  for (mod_name in names(models)) {
    # Extract the model
    mod_i <- models[[mod_name]]
    
    # Compute LOO results
    mod_loo_results_i <- loo::loo(mod_i)
    
    # Extract relevant LOO metrics
    elpd_loo_i <- mod_loo_results_i$elpd_loo
    p_loo_i <- mod_loo_results_i$p_loo
    looic_i <- mod_loo_results_i$looic
    
    # Create a data frame with metrics
    df_i <- data.frame(
      elpd_loo = elpd_loo_i,
      p_loo = p_loo_i,
      looic = looic_i,
      model = mod_name
    )
    
    est_i <- tidy(mod_i, effects = "fixed", conf.int = TRUE) %>% 
      dplyr::select(term, estimate, conf.low, conf.high) %>% 
      tidyr::pivot_wider(
        names_from = term,                # Use `term` as column names
        values_from = c(estimate, conf.low, conf.high),  # Values to pivot
        names_sep = "_"                   # Add a separator to column names
      )
    
    df_i <- cbind(df_i, est_i)
    
    # Store the data frame in the list
    loo_results_list[[mod_name]] <- df_i
  }
  
  # Combind 
  loo_results_combined <- bind_rows(loo_results_list)
  
  # Get R2 
  r2_results <- map_dfr(models, ~ as.data.frame(bayes_R2(.x)), .id = "model") %>%
    tibble::remove_rownames()
  
  # Combind R2 and loo results 
  model_fit_df <- dplyr::full_join(loo_results_combined, r2_results, by = "model") %>% 
    dplyr::select(model, everything()) %>% 
    dplyr::rename(r2 = Estimate,
                  r2_error = Est.Error,
                  r2_q2.5 = Q2.5,
                  r2_q97.5 = Q97.5) %>% 
    dplyr::mutate(id = sub("_(lm_\\d+)$", "", model),
                  model_type = sub("^.*_(lm_\\d+)$", "\\1", model))
  
  return(model_fit_df)
}

bayes_mod_predictions(): This function extracts the predicted values using fitted() from a list of models and combines them with the original data file used for the model. These are the posterior mean fitted values (i.e. the expected value of the response variable given the predictor variables and the estimated posterior distributions of the parameters) for each observation in the dataset, along with 95% credible intervals.

bayes_mod_predictions <- function(models, original_data) {

    prediction_list <- list()

    for (mod_name in names(models)) {
        # Extract mod
        mod_i <- models[[mod_name]]

        # Make mode predictions
        model_predictions_i <- as.data.frame(fitted(mod_i, summary = TRUE)) %>%
            dplyr::mutate(model = mod_name, id = sub("_(lm_\\d+)$", "", mod_name),
                model_type = sub("^.*_(lm_\\d+)$", "\\1", mod_name)) %>%
            dplyr::rename(pred_lower = Q2.5, pred_upper = Q97.5, predicted = Estimate,
                pred_error = Est.Error) %>%
            dplyr::select(model, everything())

        id_i <- model_predictions_i$id[1]

        original_data_i <- original_data %>%
            dplyr::filter(id == id_i) %>%
            dplyr::select(-id)

        model_predictions_original_i <- cbind(model_predictions_i, original_data_i)

        prediction_list[[mod_name]] <- model_predictions_original_i
    }
    predictions_df <- bind_rows(prediction_list)
    return(predictions_df)
}

calcSMR(): authored by Denis Chabot used to estimate SMR with several different methods Claireaux and Chabot (2016) ^[1]

^[1] Claireaux, G. and Chabot, D. (2016) Responses by fishes to environmental hypoxia: integration through Fry’s concept of aerobic metabolic scope. Journal of Fish Biology https://doi.org/10.1111/jfb.12833

calcSMR = function(Y, q = c(0.1, 0.15, 0.2, 0.25, 0.3), G = 1:4) {
    u = sort(Y)
    the.Mclust <- Mclust(Y, G = G)
    cl <- the.Mclust$classification
    # sometimes, the class containing SMR is not called 1 the following
    # presumes that when class 1 contains > 10% of cases, it contains SMR,
    # otherwise we take class 2
    cl2 <- as.data.frame(table(cl))
    cl2$cl <- as.numeric(levels(cl2$cl))
    valid <- cl2$Freq >= 0.1 * length(time)
    the.cl <- min(cl2$cl[valid])
    left.distr <- Y[the.Mclust$classification == the.cl]
    mlnd = the.Mclust$parameters$mean[the.cl]
    CVmlnd = sd(left.distr)/mlnd * 100
    quant = quantile(Y, q)
    low10 = mean(u[1:10])
    low10pc = mean(u[6:(5 + round(0.1 * (length(u) - 5)))])
    # remove 5 outliers, keep lowest 10% of the rest, average Herrmann & Enders
    # 2000
    return(list(mlnd = mlnd, quant = quant, low10 = low10, low10pc = low10pc, cl = cl,
        CVmlnd = CVmlnd))
}

calcO2crit(): authored by Denis Chabot used to estimate O2crit (Pcript). Claireaux and Chabot (2016) ^[1]

Note: O2 is assumed to be in percentage of dissolved oxygen (DO)

calcO2crit <- function(Data, SMR, lowestMO2 = NA, gapLimit = 4, max.nb.MO2.for.reg = 20) {
    # AUTHOR: Denis Chabot, Institut Maurice-Lamontagne, DFO, Canada first
    # version written in June 2009 last updated in January 2015
    method = "LS_reg"  # will become 'through_origin' if intercept is > 0
    if (is.na(lowestMO2))
        lowestMO2 = quantile(Data$MO2[Data$DO >= 80], p = 0.05)
    # Step 1: identify points where MO2 is proportional to DO
    geqSMR = Data$MO2 >= lowestMO2
    pivotDO = min(Data$DO[geqSMR])
    lethal = Data$DO < pivotDO
    N_under_SMR = sum(lethal)  # points available for regression?
    final_N_under_SMR = lethal  # some points may be removed at Step 4
    lastMO2reg = Data$MO2[Data$DO == pivotDO]  # last MO2 when regulating
    if (N_under_SMR > 1)
        theMod = lm(MO2 ~ DO, data = Data[lethal, ])
    # Step 2, add one or more point at or above SMR 2A, when there are fewer
    # than 3 valid points to calculate a regression
    if (N_under_SMR < 3) {
        missing = 3 - sum(lethal)
        not.lethal = Data$DO[geqSMR]
        DOlimit = max(sort(not.lethal)[1:missing])  # highest DO acceptable
        # to reach a N of 3
        addedPoints = Data$DO <= DOlimit
        lethal = lethal | addedPoints
        theMod = lm(MO2 ~ DO, data = Data[lethal, ])
    }
    # 2B, add pivotDO to the fit when Step 1 yielded 3 or more values?
    if (N_under_SMR >= 3) {
        lethalB = Data$DO <= pivotDO  # has one more value than 'lethal'
        regA = theMod
        regB = lm(MO2 ~ DO, data = Data[lethalB, ])
        large_slope_drop = (coef(regA)[2]/coef(regB)[2]) > 1.1  # arbitrary
        large_DO_gap = (max(Data$DO[lethalB]) - max(Data$DO[lethal])) > gapLimit
        tooSmallMO2 = lastMO2reg < SMR
        if (!large_slope_drop & !large_DO_gap & !tooSmallMO2)
            {
                lethal = lethalB
                theMod = regB
            }  # otherwise we do not accept the additional point
    }
    # Step 3 if the user wants to limit the number of points in the regression
    if (!is.na(max.nb.MO2.for.reg) & sum(lethal) > max.nb.MO2.for.reg) {
        Ranks = rank(Data$DO)
        lethal = Ranks <= max.nb.MO2.for.reg
        theMod = lm(MO2 ~ DO, data = Data[lethal, ])
        final_N_under_SMR = max.nb.MO2.for.reg
    }
    # Step 4
    predMO2 = as.numeric(predict(theMod, data.frame(DO = Data$DO)))
    Data$delta = (Data$MO2 - predMO2)/predMO2 * 100  # residuals set to zero
    # when below pivotDO
    Data$delta[Data$DO < pivotDO | lethal] = 0
    tol = 0  # any positive residual is unacceptable
    HighValues = Data$delta > tol
    Ranks = rank(-1 * Data$delta)
    HighMO2 = HighValues & Ranks == min(Ranks)  # keep largest residual
    if (sum(HighValues) > 0)
        {
            nblethal = sum(lethal)
            Data$W = NA
            Data$W[lethal] = 1/nblethal
            Data$W[HighMO2] = 1
            theMod = lm(MO2 ~ DO, weight = W, data = Data[lethal | HighMO2, ])
            # This new regression is always an improvement, but there can still
            # be points above the line, so we repeat
            predMO2_2 = as.numeric(predict(theMod, data.frame(DO = Data$DO)))
            Data$delta2 = (Data$MO2 - predMO2_2)/predMO2_2 * 100
            Data$delta2[Data$DO < pivotDO] = 0
            tol = Data$delta2[HighMO2]
            HighValues2 = Data$delta2 > tol
            if (sum(HighValues2) > 0)
                {
                  Ranks2 = rank(-1 * Data$delta2)
                  HighMO2_2 = HighValues2 & Ranks2 == 1  # keep the largest residual
                  nblethal = sum(lethal)
                  Data$W = NA
                  Data$W[lethal] = 1/nblethal
                  Data$W[HighMO2_2] = 1
                  theMod2 = lm(MO2 ~ DO, weight = W, data = Data[lethal | HighMO2_2,
                    ])
                  # is new slope steeper than the old one?
                  if (theMod2$coef[2] > theMod$coef[2]) {
                    theMod = theMod2
                    HighMO2 = HighMO2_2
                  }
                }  # end second search for high value
        }  # end first search for high value
    Coef = coefficients(theMod)
    # Step 5, check for positive intercept
    AboveOrigin = Coef[1] > 0
    # if it is, we use a regression that goes through the origin
    if (AboveOrigin) {
        theMod = lm(MO2 ~ DO - 1, data = Data[lethal, ])
        Coef = c(0, coefficients(theMod))  # need to add the intercept (0)
        # manually to have a pair of coefficients
        method = "through_origin"
        HighMO2 = rep(FALSE, nrow(Data))  # did not use the additional value
        # from Step 4
    }
    po2crit = as.numeric(round((SMR - Coef[1])/Coef[2], 1))
    sum_mod = summary(theMod)
    anov_mod = anova(theMod)
    O2CRIT = list(o2crit = po2crit, SMR = SMR, Nb_MO2_conforming = N_under_SMR, Nb_MO2_conf_used = final_N_under_SMR,
        High_MO2_required = sum(HighMO2) == 1, origData = Data, Method = method,
        mod = theMod, r2 = sum_mod$r.squared, P = anov_mod$"Pr(>F)", lethalPoints = which(lethal),
        AddedPoints = which(HighMO2))
}  # end function

plotO2crit(): authored by Denis Chabot, used to plot the modes used for the calcO2crit() function. Claireaux and Chabot (2016) ^[1]

Note: I added abline(h=lowestMO2, col=“pink”) so that I could visualise the lowestMO2 position

plotO2crit <- function(o2critobj, plotID = "", Xlab = "Dissolved oxygen (% sat.)",
    Ylab = "dotitalumol", smr.cex = 0.9, o2crit.cex = 0.9, plotID.cex = 1.2, Transparency = T,
    ...) {
    # AUTHOR: Denis Chabot, Institut Maurice-Lamontagne, DFO, Canada first
    # version written in June 2009 last updated 2015-02-09 for R plotting
    # devices that do not support transparency (e.g., postscript), set
    # Transparency to FALSE
    smr = o2critobj$SMR
    if (Ylab %in% c("dotitalumol", "italumol", "dotumol", "umol", "dotitalmg", "italmg",
        "dotmg", "mg")) {
        switch(Ylab, dotitalumol = {
            mo2.lab = expression(paste(italic(dot(M))[O[2]], " (", mu, "mol ", O[2],
                " ", min^-1, " ", kg^-1, ")"))
        }, italumol = {
            mo2.lab = expression(paste(italic(M)[O[2]], " (", mu, "mol ", O[2], " ",
                min^-1, " ", kg^-1, ")"))
        }, dotumol = {
            mo2.lab = expression(paste(dot(M)[O[2]], " (", mu, "mol ", O[2], " ",
                min^-1, " ", kg^-1, ")"))
        }, umol = {
            mo2.lab = expression(paste(M[O[2]], " (", mu, "mol ", O[2], " ", min^-1,
                " ", kg^-1, ")"))
        }, dotitalmg = {
            mo2.lab = expression(paste(italic(dot(M))[O[2]], " (mg ", O[2], " ",
                h^-1, " ", kg^-1, ")"))
        }, italmg = {
            mo2.lab = expression(paste(italic(M)[O[2]], " (mg ", O[2], " ", h^-1,
                " ", kg^-1, ")"))
        }, dotmg = {
            mo2.lab = expression(paste(dot(M)[O[2]], " (mg ", O[2], " ", h^-1, " ",
                kg^-1, ")"))
        }, mg = {
            mo2.lab = expression(paste(M[O[2]], " (mg ", O[2], " ", h^-1, " ", kg^-1,
                ")"))
        })
    } else mo2.lab = Ylab
    if (Transparency) {
        Col = c(rgb(0, 0, 0, 0.7), "red", "orange")
    } else {
        Col = c(grey(0.3), "red", "orange")
    }
    Data = o2critobj$origData
    lowestMO2 = quantile(Data$MO2[Data$DO >= 80], p = 0.05)  # I added this
    Data$Color = Col[1]
    Data$Color[o2critobj$lethalPoints] = Col[2]
    Data$Color[o2critobj$AddedPoints] = Col[3]
    # ordinary LS regression without added points: blue line, red symbols
    # ordinary LS regression with added points: blue line, red & orange symbols
    # regression through origin: green dotted line, red symbols
    line.color = ifelse(o2critobj$Method == "LS_reg", "blue", "darkgreen")
    line.type = ifelse(o2critobj$Method == "LS_reg", 1, 3)
    limX = c(0, max(Data$DO))
    limY = c(0, max(Data$MO2))
    plot(MO2 ~ DO, data = Data, xlim = limX, ylim = limY, col = Data$Color, xlab = Xlab,
        ylab = mo2.lab, ...)
    coord <- par("usr")
    if (plotID != "") {
        text(0, coord[4], plotID, cex = plotID.cex, adj = c(0, 1.2))
    }
    abline(h = lowestMO2, col = "pink")  # I added this
    abline(h = smr, col = "orange")
    text(coord[1], smr, "SMR", adj = c(-0.1, 1.3), cex = smr.cex)
    text(coord[1], smr, round(smr, 1), adj = c(-0.1, -0.3), cex = smr.cex)
    if (!is.na(o2critobj$o2crit)) {
        abline(o2critobj$mod, col = line.color, lty = line.type)
        segments(o2critobj$o2crit, smr, o2critobj$o2crit, coord[3], col = line.color,
            lwd = 1)
        text(x = o2critobj$o2crit, y = 0, o2critobj$o2crit, col = line.color, cex = o2crit.cex,
            adj = c(-0.1, 0.5))
    }
}  # end of function

6 📂 Directories

6.1 Input

📥 input_data_wd: Directory for the metadata

wd <- getwd()
input_data_wd <- paste0(wd, "./input-data")  # creates a variable with the name of the wd we want to use

📥 mod_data_wd: Directory for model output data estimated slopes

mod_data_wd <- paste0(wd, "./mod-data")

6.2 Output

📤 output_fig_wd: this is where we will put the figures

## [1] "Folder already exists"

📤 output_mods_wd: this is where we will put the models

## [1] "Folder already exists"

7 💿 Data

7.1 Slopes (MO₂)

💿 slope_df: We have imported the slopes extracted in LabChart during each phase of the experiment

 setwd(input_data_wd)
# 
# # Get the names of all sheets in the Excel file
sheet_names <- excel_sheets("labchart-all-dates-2025.xlsx")
all_trials_select <- c("start_date", "order", "phase", "cycle", "date", "time")
slope_list <- list()

for (sheet in sheet_names) {

  df <- read_excel("labchart-all-dates-2025.xlsx", sheet = sheet) %>% 
  dplyr::rename_with(tolower)
  
a_name <- paste0("a_", tolower(sheet))
a_df <- df %>%
  dplyr::select(starts_with('a'), all_trials_select) %>% 
  dplyr::rename(temp = a_temp) %>% 
  dplyr::mutate(across(starts_with('a'), as.numeric)) %>% 
  pivot_longer(
    cols = starts_with('a'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>%
  dplyr::mutate(respirometer_group = "a") # Add a new column with a fixed value

slope_list[[a_name]]<- a_df

b_name <- paste0("b_", tolower(sheet))
b_df <- df %>% 
  dplyr::select(starts_with('b'), all_trials_select) %>% 
  dplyr::rename(temp = b_temp) %>% 
  dplyr::mutate(across(starts_with('b'), as.numeric)) %>% 
  pivot_longer(
    cols = starts_with('b'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>% 
    dplyr::mutate(respirometer_group = "b")

slope_list[[b_name]] <- b_df

c_name <- paste0("c_", tolower(sheet))
c_df <- df %>% 
  dplyr::select(starts_with('c'), all_trials_select) %>% 
  dplyr::rename(temp = c_temp,
                i_cycle = cycle) %>% 
  dplyr::mutate(across(starts_with('c'), as.numeric)) %>%
  pivot_longer(
    cols = starts_with('c'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>% 
    dplyr::mutate(respirometer_group = "c") %>% 
  dplyr::rename(cycle = i_cycle)

slope_list[[c_name]] <- c_df

d_name <- paste0("d_", tolower(sheet))
d_df <- df %>% 
  dplyr::select(starts_with('d'), all_trials_select) %>% 
  dplyr::rename(temp = d_temp,
                i_date = date) %>% 
  dplyr::mutate(across(starts_with('d'), as.numeric)) %>%
  pivot_longer(
    cols = starts_with('d'), # Select all columns to pivot
    names_to = c("chamber_id", ".value"), # Separate column names into 'id' and other variables
    names_sep = "_"
  ) %>% 
    dplyr::mutate(respirometer_group = "d") %>% 
  dplyr::rename(date = i_date)

slope_list[[d_name]] <- d_df
}


slope_df <- bind_rows(slope_list) %>% 
  dplyr::mutate(resp_cat_date = paste0(respirometer_group, "_", start_date),
                chamber_n = str_extract(chamber_id, "\\d+"),
                id_prox = paste0(resp_cat_date, "_", chamber_n),
                time_hms = as_hms(time*3600),
                date_chr = format(date, "%d/%m/%Y")
                )

7.2 Fish metadata

💿 metadata: This is the meta data for each chamber

Note: We are also adding volume based on chamber type.

setwd(input_data_wd)

metadata <- read_excel("fish-meta.xlsx", na = "NA") %>%
    dplyr::mutate(id_split = id) %>%
    tidyr::separate(id_split, into = c("respirometer_group", "salinity_group", "start_date",
        "chamber"), sep = "_") %>%
    dplyr::mutate(volume = dplyr::case_when(chamber_type == "L" ~ 0.3, chamber_type ==
        "M_M" ~ 0.105, chamber_type == "M_NM" ~ 0.11, chamber_type == "S" ~ 0.058,
        chamber_type == "SM" ~ 0.075, chamber_type == "D3" ~ 0.055, TRUE ~ NA), id_prox = paste0(respirometer_group,
        "_", start_date, "_", chamber), rel_size = mass/volume)

7.2.1 Urbina, Glover, and Forster (2012)

💿 urbina_et_al_2012: This is the mean level data extracted from Urbina, Glover, and Forster (2012)^[2] Figure 1a. We used the metaDigitise package in R to extract the data ^[3].

setwd(input_data_wd)

urbina_et_al_2012 <- read.csv("urbina-et-al-2012-fig1a.csv")

7.2.2 Pcrit viusal inspection

💿 pcrit_check: This data frame is a list of all fish that had higher order polynomials (2^nd and 3^rd order; see ‘Incremental regression analyses’ for details) and subsequently had Pcrit modelled (see ‘Pcrit model’ for details). In this data we have checked each Pcrit model and confirmed visual if a Pcrit is present.

setwd(input_data_wd)

pcrit_check <- read_excel("pcrit-check.xlsx", na = "NA")

7.2.3 Combinding metadata

Adding the meta data to the slopes data frame

slope_df_2 <- slope_df %>%
    dplyr::select(-start_date, -respirometer_group) %>%
    left_join(metadata, by = "id_prox") %>%
    dplyr::mutate(light_dark = if_else(time_hms >= as.hms("07:00:00") & time_hms <
        as.hms("19:00:00"), "light", "dark")) %>%
    dplyr::arrange(id)

8 🧼 Data tidy

8.1 Numbers

We have 64 fish with MO₂ data

n <- slope_df_2 %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::distinct(id) %>%
    nrow(.)

paste0("n = ", n)

## [1] "n = 64"

With 48 from the 0 ppt and 48 from 9 ppt groups

slope_df_2 %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::group_by(salinity_group) %>%
    dplyr::reframe(`n total` = length(unique(id))) %>%
    gt() %>%
    cols_label(salinity_group = "Salinity group") %>%
    cols_align(align = "center", columns = everything())

Salinity group	n total
0	33
9	31

8.2 Fish size

Here we calculate the mean length and size of fish used in the experiment.

mass_length <- slope_df_2 %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::group_by(id) %>%
    dplyr::sample_n(1) %>%
    dplyr::ungroup() %>%
    dplyr::reframe(x_mass = round(mean(mass, na.rm = TRUE), 3), min_mass = round(min(mass,
        na.rm = TRUE), 3), max_mass = round(max(mass, na.rm = TRUE), 3), x_length = round(mean(length,
        na.rm = TRUE), 2), min_length = round(min(length, na.rm = TRUE), 2), max_length = round(max(length,
        na.rm = TRUE), 2))

mass_mean <- mass_length %>%
    pull(x_mass)

mass_min <- mass_length %>%
    pull(min_mass)

mass_max <- mass_length %>%
    pull(max_mass)

length_mean <- mass_length %>%
    pull(x_length)

length_min <- mass_length %>%
    pull(min_length)

length_max <- mass_length %>%
    pull(max_length)

paste0("The mean mass of fish was ", mass_mean, " g (range: ", mass_min, "-", mass_max,
    ")", ", and the mean length was ", length_mean, " mm (range: ", length_min, "-",
    length_max, ")")

## [1] "The mean mass of fish was 0.532 g (range: 0.21-1.6), and the mean length was 50.41 mm (range: 40-70)"

8.3 Filtering trials

We will remove 6 trials which had errors. These are as follows:

• a_0_25nov_3 fish died during trial
• b_0_26nov_4 flat lined early
• c_0_22nov_2 chamber was opened early
• c_9_26nov_2 chamber was opened early
• c_9_26nov_4 chamber was opened early
• d_9_27nov_3 sensor was jumpy and end points were hard to confidently ID visually

remove_trial_error <- c("a_0_25nov_3", "b_0_26nov_4", "c_0_22nov_2", "c_9_26nov_2",
    "c_9_26nov_4", "d_9_27nov_3")

slope_df_filter <- slope_df_2 %>%
    dplyr::filter(!(id %in% remove_trial_error))

We now have 58 fish with MO2 data

n <- slope_df_filter %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::distinct(id) %>%
    nrow(.)

paste0("n = ", n)

## [1] "n = 58"

With 30 in the 0 ppt group and 28 in the 9 ppt group

slope_df_filter %>%
    dplyr::filter(chamber_condition == "fish") %>%
    dplyr::group_by(salinity_group) %>%
    dplyr::reframe(`n total` = length(unique(id))) %>%
    gt() %>%
    cols_label(salinity_group = "Salinity group") %>%
    cols_align(align = "center", columns = everything())

Salinity group	n total
0	30
9	28

8.4 Filtering MO₂ estimates

Here we apply the following filters to the MO₂ data:

• Remove the first 5 SMR cycles (burn in)
• Remove all positive raw slopes
• Remove all MO₂ calculated using less then 60 data points (5 min)
• Remove all MO₂ calculated if O₂ increases in a closed phase (i.e. trial has ended)

Check positive values for MO₂ before removing.

slope_tidy_remove_flush <- slope_df_filter %>%
  dplyr::filter(phase != "smr", chamber_condition == "fish") %>% 
  dplyr::group_by(id) %>%
  dplyr::arrange(id, order) %>%  # Ensure the data is sorted within each group
  dplyr::mutate(o2_diff = if_else(row_number() == 1, 0, o2 - lag(o2)), # Calculate the difference in 'o2'
                o2_diff_cumsum = cumsum(o2_diff > 1)) %>%  # Checks first occurrence and sums 
  dplyr::filter(o2_diff_cumsum == 0) %>%  # Keep rows until the first jump > 1
  dplyr::ungroup() %>%
  dplyr::select(-o2_diff, -o2_diff_cumsum)

postive_slopes <- slope_tidy_remove_flush %>%
  dplyr::filter(chamber_condition == "fish" & mo2corr > 0)

list_postive_all <- slope_df_filter %>% 
  dplyr::filter(chamber_condition == "fish" & mo2corr > 0 & phase == "smr") %>% 
  dplyr::bind_rows(., postive_slopes) %>% 
  dplyr::distinct(id) %>% 
  dplyr::pull(id)

print(paste0("There are ", length(list_postive_all), " fish with postive slopes. These fish are: ", paste(list_postive_all, collapse = ", "), "."))

## [1] "There are 22 fish with postive slopes. These fish are: a_9_22nov_1, d_0_21nov_3, a_0_24nov_1, a_0_25nov_1, b_0_24nov_1, b_0_24nov_3, b_0_25nov_1, b_0_25nov_4, b_9_21nov_1, b_9_21nov_2, b_9_21nov_3, b_9_22nov_1, b_9_22nov_2, b_9_22nov_3, c_0_21nov_1, c_0_21nov_2, c_0_22nov_4, c_9_24nov_3, c_9_25nov_3, d_0_22nov_2, d_9_25nov_2, d_9_25nov_3."

Filtering the MO₂ data

cycle_burn <- 0:4

slope_df_filter_1 <- slope_df_filter %>%
  dplyr::filter(!(cycle %in% cycle_burn) & 
                  mo2corr < 0 & 
                  n > 60 &
                  chamber_condition == "fish"
                )
  
# Filter out the end flush
slope_tidy_closed <- slope_df_filter_1 %>%
  dplyr::filter(phase != "smr") %>% 
  dplyr::group_by(id) %>%
  dplyr::arrange(id, order) %>%  # Ensure the data is sorted within each group
  dplyr::mutate(o2_diff = if_else(row_number() == 1, 0, o2 - lag(o2)), # Calculate the difference in 'o2'
                o2_diff_cumsum = cumsum(o2_diff > 1)) %>%  # Checks first occurrence and sums 
  dplyr::filter(o2_diff_cumsum == 0) %>%  # Keep rows until the first jump > 1
  dplyr::ungroup() %>%
  dplyr::select(-o2_diff, -o2_diff_cumsum)

slope_tidy_smr <- slope_df_filter_1 %>% 
  dplyr::filter(phase == "smr")

slope_df_filter_2 <- rbind(slope_tidy_smr, slope_tidy_closed) %>% 
  dplyr::arrange(id, order)

8.5 Calculating SMR

We will calculate SMR using calcSMR function by Chabot, Steffensen and Farrell (2016)^[1]. Specifically, we use mean of the lowest normal distribution (MLND) where CVmlnd < 5.4 and the mean of the lower 20% quantile (q0.2) were CVmlnd > 5.4. If CVmlnd is not calculated we have used q0.2.

labchart_chabot_smr <- slope_df_filter_2 %>%
    dplyr::filter(phase == "smr")

# Extract distinct IDs
ids <- labchart_chabot_smr %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

# Initialise an empty list to store SMR data
smr_list <- list()

# Process each ID
for (id_i in ids) {
    tryCatch({
        # Filter the data for the current ID
        df_i <- labchart_chabot_smr %>%
            dplyr::filter(id == id_i) %>%
            dplyr::mutate(abs_mo2corr = abs(mo2corr))

        # Calculate SMR results
        calcSMR_results <- calcSMR(df_i$abs_mo2corr)
        CVmlnd_i <- calcSMR_results$CVmlnd
        quant_i <- calcSMR_results$quant %>%
            as_tibble()
        quant_20per_i <- quant_i$value[3]
        mlnd_i <- calcSMR_results$mlnd
        smr_value <- if_else(CVmlnd_i < 5.4, mlnd_i, quant_20per_i)
        smr_type <- if_else(CVmlnd_i < 5.4, "mlnd", "quant_20per")
        smr_value <- if_else(is.na(smr_value), quant_20per_i, smr_value)
        smr_type <- if_else(is.na(smr_type), "quant_20per", smr_type)

        # Create a data frame for the current ID
        smr_df <- tibble(id = id_i, smr = smr_value, smr_est = smr_type)

    }, error = function(e) {
        # Handle errors by assigning NA values
        smr_df <- tibble(id = id_i, smr = NA, smr_est = NA)
    })

    # Append to the list
    smr_list[[id_i]] <- smr_df
}

# Combine all individual SMR data frames into one
smr_df <- bind_rows(smr_list) %>%
    dplyr::rename(smr_chabot = smr, smr_chabot_method = smr_est)

slope_df_filter_3 <- slope_df_filter_2 %>%
    dplyr::left_join(., smr_df, by = "id")

8.6 Transforming MO₂

Here we are transforming the MO₂ units. The resulting values are as follows:

• MO2 is absolute value of the background and leak corrected MO₂ slope from Labchart (mo2corr) times the net volume of the chamber (volume - fish mass), × 60, × 60, to achieve MO₂ as mg^-1 O₂ h^-1

• MO2_g is MO2 divided by fish mass to achieve MO₂ as mg^-1 O₂ g^-1 h^-1 (i.e. mass standardised)

• SMR absolute value of the SMR estimates using methods described by Chabot, Steffensen and Farrell (2016)^[1] times the net volume of the chamber (volume - fish mass), × 60, × 60, to achieve SMR as mg^-1 O₂ g^-1 h^-1)

• SMR_g is SMR divided by fish mass to achieve SMR as mg^-1 O₂ g^-1 h^-1 (i.e. mass standardised)

• DO is dissolved oxygen percentage calculated from O₂ values (mg^-1 L^-1) using the recorded temperature, salinity, and a constant atmospheric pressure (Pa; 1013.25)

• o2_kpa is the O₂ concentration in kilopascal (kpa). This is used to make a comparative figure only.

# Combine back into one data frame
slope_tidy <- slope_df_filter_3 %>% 
    dplyr::mutate(DO = conv_o2(
                   o2 = o2,
                   from = "mg_per_l",
                   to = "percent_a.s.",
                   temp = temp, #C
                   sal = measured_salinity,
                   atm_pres = 1013.25),
                  o2_kpa = conv_o2(
                   o2 = o2,
                   from = "mg_per_l",
                   to = "kPa",
                   temp = temp, #C
                   sal = measured_salinity,
                   atm_pres = 1013.25),
                  net_volume = volume - mass/1000,
                  MO2 = abs(mo2corr)*net_volume*60*60,
                  MO2_g = MO2/mass,
                  SMR = abs(smr_chabot)*net_volume*60*60,
                  SMR_g = SMR/mass
                  )

9 📊 Visualising data

9.1 O₂ vs MO₂

9.1.1 Figure S1

This interactive was used to identify any outliers, or potential errors.

lm_lines <- slope_tidy %>%
  group_by(id) %>%
  summarise(
    DO_seq = list(seq(min(DO), max(DO), length.out = 100)),  # Generate a sequence of DO values
    MO2_pred = list(predict(lm(MO2_g ~ DO, data = cur_data()), newdata = data.frame(DO = seq(min(DO), max(DO), length.out = 100))))
  ) %>%
  unnest(c(DO_seq, MO2_pred))  # Expand lists into rows

# Create scatter plot with markers for each fish
p <- plot_ly(
  data = slope_tidy,
  x = ~DO,
  y = ~MO2_g,
  type = "scatter",
  mode = "markers",
  color = ~id,  # Colour points by fish ID
  marker = list(opacity = 0.6),
  name = ~id
)

# Add regression lines for each fish
p <- p %>%
  add_trace(
    data = lm_lines,
    x = ~DO_seq,
    y = ~MO2_pred,
    type = "scatter",
    mode = "lines",
    color = ~id,  # Ensure each line matches its corresponding fish
    line = list(width = 1, dash = "solid"),
    showlegend = FALSE  # Avoid cluttering the legend
  )

# Final layout
p <- p %>%
  layout(
    title = "MO<sub>2</sub> vs Dissolved Oxygen with individual linear regressions",
    xaxis = list(title = "Dissolved Oxygen (%)"),
    yaxis = list(title = "MO<sub>2</sub> (mg<sup>-1</sup> O<sub>2</sub> g<sup>-1</sup> h<sup>-1</sup>)"),
    showlegend = FALSE
  )

# Display plot
p

Figure S1: interactive plot of metabolic rate measurements (MO₂; mg O₂ g^-1h^-1) by dissolved oxygen percentage (DO) for all fish, including all estimates during the SMR phase (i.e. intermittent phase). Individual linear regression were fitted for visual reference, and do not represent the best fitting regression.

9.1.2 Figure S2

Looking at the difference responses in the two salinity groups.

salinity_summary <- slope_tidy %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(n = length(unique(id)))

slope_tidy %>% 
  ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
  geom_point(show.legend = FALSE) +
  geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
  geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
  geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs"), se = TRUE, colour = "red", linetype = "dashed") +
  theme_clean() +
  facet_wrap(~salinity_group) +
  labs(
    subtitle = "mo2 vs o2 by salinity treatment",
    x = "Dissolved oxygen percentage (DO)",
    y = "MO2 (O2 mg/g/h)"
  ) +
  geom_text(data = salinity_summary,
            aes(x = -Inf, y = Inf, label = paste0("italic(n) == ", n)),
            hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, parse = TRUE)

Figure S2: Metabolic rate measurements (MO₂; mg O₂ g^-1h^-1) by dissolved oxygen percentage (DO) for fish from the two salinity treatments, including all estimates during the SMR phase (i.e. intermittent phase). Individual linear regression were fitted for visual reference, and do not represent the best fitting regression. The solid red and dashed red line represents the global trend in each of the treatments, both a linear (solid red) and non-linear (dashed red; Generalized Additive Model fitted with geom_smooth).

9.1.3 Figure S3

A plot to look at the different chamber types

chamber_summary <- slope_tidy %>% 
  dplyr::group_by(chamber_type) %>% 
  dplyr::reframe(n = length(unique(id)))

slope_tidy %>% 
  ggplot(aes(y = MO2_g, x = DO, colour = id)) + # Default aesthetics
  geom_point(show.legend = FALSE) +
  geom_smooth(aes(group = id), method = "lm", se = FALSE, colour = scales::alpha("black", 0.5)) + # Transparent black lines
  geom_smooth(method = "lm", se = TRUE, colour = "red") + # Overall smooth line
  geom_smooth(se = TRUE, colour = "red", linetype = "dashed") +
  theme_clean() +
  facet_wrap(~chamber_type, scale = "free") +
  labs(
    subtitle = "mo2 vs o2 by chamber type",
    x = "Dissolved oxygen percentage (DO)",
    y = "MO2 (mg O2 g/h)"
  ) +
  geom_text(data = chamber_summary,
            aes(x = -Inf, y = Inf, label = paste0("italic(n) == ", n)),
            hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, parse = TRUE)

Figure S3: Metabolic rate measurements (MO₂; mg O₂ g^-1h^-1) by dissolved oxygen percentage (DO) for fish tested in the 4 different chamber types, including all estimates during the SMR phase (i.e. intermittent phase). Individual linear regression were fitted for visual reference, and do not represent the best fitting regression. The solid red and dashed red line represents the global trend in each of the treatments, both a linear (solid red) and non-linear (dashed red; Generalized Additive Model fitted with geom_smooth).

9.1.4 Figure S4

Here, we have recreated a similar plot to that presented in Urbina, Glover, and Forster (2012)^[1] and have extracted the mean level data from Figure 1a using the metaDigitise package in R^[3]. This data is called urbina_et_al_2012. This allows us to compare the differences in Mo2 and the relationship between Mo2 and O₂.

First making a binned data frame to match Urbina, Glover, and Forster (2012) as closely as possible.

min_o2_kpa <- min(slope_tidy$o2_kpa, na.rm = TRUE)
max_o2_kpa <- max(slope_tidy$o2_kpa, na.rm = TRUE)

time_bin_df <- slope_tidy %>%
  mutate(o2_group = cut(o2_kpa, 
                        breaks = seq(min_o2_kpa, max_o2_kpa, length.out = 13), # 11 intervals, so 12 breakpoints
                        labels = paste0("Group ", 1:12), 
                        include.lowest = TRUE)) %>% 
  dplyr::group_by(o2_group) %>% 
  dplyr::reframe(mean_MO2_g = mean(MO2_g)*31.25,
                 mean_o2_kpa = mean(o2_kpa),
                 n = length(MO2_g)*31.25,
                 MO2_g_sd = sd(MO2_g)*31.25,
                 o2_kpa_sd = sd(o2_kpa))

Now the plot with our mean data and the mean data from Urbina, Glover, and Forster (2012)^[1]

n <- slope_tidy %>%
    dplyr::distinct(id) %>%
    nrow(.)

# Existing plot
p <- time_bin_df %>%
    ggplot() + geom_point(data = slope_tidy, aes(y = MO2_g * 31.25, x = o2_kpa),
    size = 2, color = "grey", alpha = 0.3) + geom_point(data = time_bin_df, aes(y = mean_MO2_g,
    x = mean_o2_kpa), size = 3, colour = "#0E4C92", show.legend = FALSE) + geom_errorbar(data = time_bin_df,
    aes(y = mean_MO2_g, x = mean_o2_kpa, ymin = mean_MO2_g - MO2_g_sd, ymax = mean_MO2_g +
        MO2_g_sd), width = 0.15, colour = "#0E4C92") + geom_errorbarh(data = time_bin_df,
    aes(y = mean_MO2_g, x = mean_o2_kpa, xmin = mean_o2_kpa - o2_kpa_sd, xmax = mean_o2_kpa +
        o2_kpa_sd), height = 0.4, colour = "#0E4C92") + geom_point(data = urbina_et_al_2012,
    aes(x = o2_mean, y = mo2_mean), size = 3, shape = 21, fill = "#D21F3C", color = "#D21F3C",
    stroke = 1) + geom_errorbar(data = urbina_et_al_2012, aes(x = o2_mean, ymin = mo2_mean -
    mo2_sd, ymax = mo2_mean + mo2_sd), width = 0.2, colour = "#D21F3C") + geom_errorbarh(data = urbina_et_al_2012,
    aes(y = mo2_mean, xmin = o2_mean - o2_sd, xmax = o2_mean + o2_sd), height = 0.4,
    colour = "#D21F3C") + annotate("text", x = 0, y = 16, label = bquote(atop("Current data (blue), " *
    italic(n) * " = " * .(n), "Urbina data (red), " * italic(n) * " = 67")), hjust = 0,
    vjust = 1, size = 4) + theme_clean() + labs(subtitle = "", x = "PO2 (kPa)", y = "MO2 (umol O2 g/h)") +
    scale_y_continuous(limits = c(0, 16), breaks = seq(0, 16, by = 2)) + scale_x_continuous(limits = c(0,
    22), breaks = seq(0, 22, by = 2))

p

Figure S4: Mean and standard error of metabolic rate (MO₂ μmol O₂ g^-1 h^-1) and oxygen concentration (PO₂ kPa) using 12 evenly spaced bins over the O₂ range of observed data (blue dots). Compared against the mean and standard error reported in Urbina (2012)^[1] (red dots). The grey dots are the raw observed data form the present study.

9.2 Routine MO₂

Making an SMR phase only data frame

slope_tidy_smr <- slope_tidy %>%
    dplyr::filter(phase == "smr")

9.2.1 Figure S5

Routine MO₂ by salinity

mean_mo2_salinity <- slope_tidy_smr %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(mean_mo2 = mean(MO2, na.rm = TRUE))

fig_i <- ggplot() +
    geom_violin(data = slope_tidy_smr, aes(x = salinity_group, y = MO2, fill = salinity_group), color = NA, alpha = 0.3) +
  geom_jitter(data = slope_tidy_smr, aes(x = salinity_group, y = MO2, fill = salinity_group),
                       shape = 21, size = 2, color = "black", alpha = 0.2) +
  geom_boxplot(data = slope_tidy_smr, aes(x = salinity_group, y = MO2, fill = salinity_group),
                        size = 1, alpha = 0.5, outlier.shape = NA, width = 0.3) +
  geom_point(data = mean_mo2_salinity, 
                aes(x = salinity_group, y = mean_mo2, fill = salinity_group), 
                size = 3, alpha = 0.8, colour = "black", stroke = 2) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme_clean() +
  theme(legend.position = "none") +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "Routine MO2 (mg O2 g/h)"
  )

fig_i

Figure S5: Plot of all MO₂ measures during SMR phase by salinity treatment. The small points are the raw observed values, the shaded area behind the points is the a kernel density of the observed data, the box plot shows the median and interquartile range (IQR), and the large point shows the mean.

9.3 SMR

9.3.1 Figure S6

Here’s the same plot but for only the SMR, as estimated with calcSMR function by Chabot, Steffensen and Farrell (2016)^[1]. Specifically, we use mean of the lowest normal distribution (MLND) where CVmlnd < 5.4 and the mean of the lower 20% quantile (q0.2) were CVmlnd > 5.4. If CVmlnd is not calculated we have used q0.2.

smr_only <- slope_tidy_smr %>% 
  dplyr::group_by(id) %>% 
  dplyr::slice(1)

mean_smr_salinity <- smr_only %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(mean_smr = mean(SMR, na.rm = TRUE))

fig_i <- ggplot() +
    geom_violin(data = smr_only, aes(x = salinity_group, y = SMR, fill = salinity_group), color = NA, alpha = 0.3) +
  geom_jitter(data = smr_only, aes(x = salinity_group, y = SMR, fill = salinity_group),
                       shape = 21, size = 2, color = "black", alpha = 0.2) +
  geom_boxplot(data = smr_only, aes(x = salinity_group, y = SMR, fill = salinity_group),
                        size = 1, alpha = 0.5, outlier.shape = NA, width = 0.3) +
  geom_point(data = mean_smr_salinity, 
                aes(x = salinity_group, y = mean_smr, fill = salinity_group), 
                size = 3, alpha = 0.8, colour = "black", stroke = 2) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme_clean() +
  theme(legend.position = "none") +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "MO2 (mg O2 g/h)"
  )

fig_i

Figure S6: The standard metabolic rate (SMR) by salinity treatment. The small points are the raw observed values, the shaded area behind the points is the a kernel density of the observed data, the box plot shows the median and interquartile range (IQR), and the large point shows the mean.

9.4 Individual O₂, MO₂, and SMR

Here we will plot the individual relationship between O₂, MO₂, and SMR for all fish

Create output directory if needed

9.4.1 Figure S7

Loop through each fish ID to create a plot, save these to a single PDF file

library(ggplot2)
library(dplyr)

# Define consistent colors for all phases
phase_colors <- c(
  "50c"  = "#e56b6f",
  "75c"  = "#b56576",
  "100c" = "#6d597a",
  "smr"  = "#355070"
)

ids <- slope_tidy %>% 
  dplyr::distinct(id) %>% 
  pull(id) %>% 
  as.list()

MO2_plot_list <- list()

# 1) Open the PDF device once
pdf(
  file   = file.path(output_fig_slopes_wd, "combined_slopes.pdf"), 
  width  = 8, 
  height = 6
)

# 2) Loop over IDs and create each plot
for (id_i in ids) {
  
  smr <- slope_tidy %>% 
    dplyr::filter(id == id_i) %>% 
    dplyr::slice(1) %>% 
    dplyr::pull(SMR)
  
  plot <- slope_tidy %>% 
    dplyr::filter(id == id_i) %>% 
    ggplot(aes(x = o2, y = MO2)) +
    geom_hline(yintercept = smr, linetype = "dashed", color = "darkred") +
    geom_point(aes(colour = phase)) +
    scale_color_manual(values = phase_colors, drop = FALSE) +  # Ensures consistent colours
    theme_clean() +
    labs(
      subtitle = paste0(id_i, " slopes"),
      x = "Mean o2 (mg_per_l)",
      y = "abs(mo2) (mg_per_l)"
    )
  
  # Instead of saving each plot separately, just print it
  print(plot)
  
  MO2_plot_list[[id_i]] <- plot
}

# 3) Close the PDF device *after* the loop
dev.off()

## png 
##   2

for (p in MO2_plot_list) {
  print(p)
}

Figure S7: Metabolic rate measurements (MO₂; mg O₂ g^-1h^-1 by O₂) at each of the four experimental phase (the over night SMR intermittent-flow respirometry phase, closed phases at 50%, 75% or 100% O₂). The estimated SMR value is represented as a dashed red line.

10 🧮 Analysis

10.1 Routine MO2

10.1.1 Scaling predictors

Here we scale our predictors for the model

scale_list <- c("temp", "order", "mass")

slope_tidy_smr <- slope_tidy_smr %>%
    dplyr::mutate(across(all_of(scale_list), ~scale(.x, center = TRUE, scale = FALSE),
        .names = "{.col}_z"), light_dark_c = if_else(light_dark == "light", 0.5,
        -0.5))

10.1.2 Model structure

Here we will use a Bayesian Generalised Linear Mixed Model (GLMM) with a Gamma distribution and a log link, where the shape parameter (K) is also modelled as a function of predictors. This models MO₂ by salinity during the SMR phase to see if the fish held at different salinities have different SMRs. We have also added a few scaled predictors, that may help describe variation in the data, such as mass (g; 0.21–1.6) temperature (°C; 13.841–14.277), measurement order (1–28), and light/dark cycle (light or dark; light between 07:00:00 and 19:00:00), we also include a random effect for fish id to account for multiple MO₂ measures on each fish (1 - 58). We allowed the the shape parameter (K) to vary as a function of some of the predictors (e.g. salinity_group, order_z) to improve fit.

mo2_gamma_bf <- bf(MO2 ~ temp_z + order_z + light_dark_c + mass_z + salinity_group +
    (1 | id), shape ~ salinity_group + order_z, family = Gamma(link = "log"))

10.1.3 Prior selection

These are the default priors. We will use these.

default_prior <- get_prior(mo2_gamma_bf, data = slope_tidy_smr, family = Gamma(link = "log"))
default_prior %>%
    gt()

prior	class	coef	group	resp	dpar	nlpar	lb	ub	source
	b								default
	b	light_dark_c							default
	b	mass_z							default
	b	order_z							default
	b	salinity_group9							default
	b	temp_z							default
student_t(3, -2.7, 2.5)	Intercept								default
student_t(3, 0, 2.5)	sd						0		default
	sd		id						default
	sd	Intercept	id						default
	b				shape				default
	b	order_z			shape				default
	b	salinity_group9			shape				default
student_t(3, 0, 2.5)	Intercept				shape				default

10.1.4 Run model

Here we run the model, I have hashed this out because I have saved the model for quick reloading.

⏭️ Skip this step if you have downloaded the saved model ‘mo2_mod_gamma.rds’

Here we reload the model

setwd(output_mods_wd)

mo2_mod_gamma <- readRDS(file = "mo2_mod_gamma.rds")

10.1.5 Model diagnostics

Checking model convergence

color_scheme_set("red")
plot(mo2_mod_gamma, ask = F)

Checking rhat are equal to one

summary(mo2_mod_gamma)

##  Family: gamma 
##   Links: mu = log; shape = log 
## Formula: MO2 ~ temp_z + order_z + light_dark_c + mass_z + salinity_group + (1 | id) 
##          shape ~ salinity_group + order_z
##    Data: slope_tidy_smr (Number of observations: 893) 
##   Draws: 4 chains, each with iter = 8000; warmup = 1000; thin = 2;
##          total post-warmup draws = 14000
## 
## Multilevel Hyperparameters:
## ~id (Number of levels: 58) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.34      0.03     0.28     0.41 1.00     3864     5351
## 
## Regression Coefficients:
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                -2.55      0.06    -2.68    -2.43 1.00     2450
## shape_Intercept           2.57      0.07     2.43     2.70 1.00    12540
## temp_z                   -0.36      0.22    -0.80     0.08 1.00     6496
## order_z                   0.00      0.00    -0.00     0.01 1.00    11555
## light_dark_c              0.14      0.02     0.09     0.19 1.00    13253
## mass_z                    1.24      0.17     0.90     1.58 1.00     3339
## salinity_group9          -0.12      0.09    -0.30     0.07 1.00     2440
## shape_salinity_group9     0.11      0.10    -0.09     0.30 1.00    12382
## shape_order_z            -0.04      0.01    -0.06    -0.02 1.00    12079
##                       Tail_ESS
## Intercept                 4218
## shape_Intercept          12262
## temp_z                    9064
## order_z                  12406
## light_dark_c             12316
## mass_z                    5807
## salinity_group9           4288
## shape_salinity_group9    12635
## shape_order_z            12068
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Using leave one out (loo) measure of fit, the model appears to preform well, all Pareto k estimates are good (k < 0.7)

loo(mo2_mod_gamma)

## 
## Computed from 14000 by 893 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   2251.6 39.0
## p_loo        71.5  7.0
## looic     -4503.3 77.9
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.0]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.

Model predictions generally align with the observed data

color_scheme_set("red")
p <- pp_check(mo2_mod_gamma, type = "dens_overlay")
p

10.1.6 📈 Results

We did not see a meaningful difference between the routine metabolic rate for fish from the two salinity treatments.

10.1.6.1 Table S1

Table S1: Fixed effect Estimates (β) and 95% Credible Intervals (95% CI)

model_est <- fixef(mo2_mod_gamma, probs = c(0.025, 0.975)) %>%
    as.data.frame() %>%
    tibble::rownames_to_column(var = "Predictor") %>%
    dplyr::mutate(β = round(Estimate, 3), Q2.5 = round(Q2.5, 3), Q97.5 = round(Q97.5,
        3), `95% CI` = paste0("[", Q2.5, ", ", Q97.5, "]"))

model_est %>%
    dplyr::select(Predictor, "β", "95% CI") %>%
    gt()

Predictor	β	95% CI
Intercept	-2.554	[-2.68, -2.427]
shape_Intercept	2.570	[2.434, 2.704]
temp_z	-0.360	[-0.804, 0.077]
order_z	0.002	[-0.002, 0.005]
light_dark_c	0.137	[0.089, 0.186]
mass_z	1.241	[0.903, 1.585]
salinity_group9	-0.115	[-0.301, 0.067]
shape_salinity_group9	0.107	[-0.087, 0.304]
shape_order_z	-0.040	[-0.059, -0.022]

Looking at the marginal mean difference between salinity groups

em_results <- emmeans(mo2_mod_gamma, ~salinity_group)
contrast_results <- contrast(em_results, method = "pairwise")
em_results_df <- em_results %>%
    tidy() %>%
    mutate(across(where(is.numeric), ~exp(.)))
contrast_results_df <- contrast_results %>%
    tidy() %>%
    mutate(across(where(is.numeric), ~exp(.)))


em_results_df %>%
    gt()

salinity_group	estimate	lower.HPD	upper.HPD
0	0.07775342	0.06851859	0.08822464
9	0.06929211	0.06056439	0.07873088

Pulling the emmeans draws for our plot

emmeans_draws <- mo2_mod_gamma %>%
    emmeans(~salinity_group) %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value), salinity_group = as.character(salinity_group))

emmeans_contrast_draws <- mo2_mod_gamma %>%
    emmeans(~salinity_group) %>%
    contrast(method = "pairwise") %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value))

10.1.6.2 Figure S8

NOTE: This plot is in the main text of the manuscript as Figure 1a

mean_mo2_salinity <- slope_tidy_smr %>% 
  dplyr::group_by(salinity_group) %>% 
  dplyr::reframe(mean_mo2 = mean(MO2, na.rm = TRUE))


fig_1 <- ggplot() +
  geom_violin(data = slope_tidy_smr,
              aes(x = salinity_group, y = MO2, fill = salinity_group),
              color = NA, alpha = 0.2) +
  geom_jitter(data = slope_tidy_smr,
              aes(x = salinity_group, y = MO2, fill = salinity_group),
              shape = 21, width = 0.3, size = 1, color = "black", alpha = 0.1) +
    geom_point(data = mean_mo2_salinity,
             aes(x = salinity_group, y = mean_mo2, fill = salinity_group),
             size = 4, alpha = 1, stroke = 2, color = "black", shape = 21,
              position = position_nudge(x = 0.05)) +
  stat_pointinterval(data = emmeans_draws, 
                     aes(x = salinity_group, y = .value),
                     color = "black", fill = "grey", point_interval = "mean_qi", .width = 0.95, shape = 21,  stroke = 2, point_size = 4, alpha = 1,
                     position = position_nudge(x = -0.05)) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme_clean() +
  theme(legend.position = "none") +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "Routine MO2 (mg O2 g/h)"
  )

fig_1

Figure 1a: Routine metabolic rate (i.e. MO₂ (mg^-1 O₂ h^-1) measured during SMR meassurments) plotted by salinity treatment. The small transparent points are the observed values, the shaded area behind the points is the a kernel density of the observed data, the large coloured point (to the right) is the observed mean, the large grey point with error bars (to the left) is the model estimated marginal mean (eemean) 95% Credible Intervals (95% CI).

10.2 SMR

10.2.1 Formating and scaling data

Here we are filtering the data frame to have only measure per fish for the SMR estimate

scale_list <- c("temp_mean", "mass", "cycles")

smr <- slope_tidy_smr %>%
    dplyr::group_by(id) %>%
    dplyr::reframe(temp_mean = mean(temp), mass = mass[1], SMR = SMR[1], salinity_group = salinity_group[1],
        cycles = length(order)) %>%
    dplyr::mutate(across(all_of(scale_list), ~scale(.x, center = TRUE, scale = FALSE),
        .names = "{.col}_z"))

10.2.2 Model structure

Here we will use a Bayesian Generalised Linear Mixed Model (GLMM) with a Gamma distribution and a log link Gamma(link = "log"), where the shape parameter (K) is also modelled as a function of the salinity group, shape ~ salinity_group. This models MO2 by salinity during the SMR phase to see if the fish held at different salinities have different SMRs. We have also added a few scaled predictors, that may help describe variation in the data, such as mass (g; 0.21–1.6) temperature (°C; 13.841–14.277), measurement order (1–28), and light/dark cycle (light or dark; light between 07:00:00 and 19:00:00), we also include a random effect for fish id to account for multiple MO2 measures on each fish (1 - 58). We allowed the the shape parameter (K) to vary as a function of some of the predictors (e.g. salinity_group, order_z) to improve fit.

smr_gamma_bf <- bf(SMR ~ temp_mean_z + cycles_z + mass_z + salinity_group, shape ~
    salinity_group, family = Gamma(link = "log"))

10.2.3 Prior selection

These are the default priors for the model. We will use these.

priors_default <- get_prior(smr_gamma_bf, data = smr, family = Gamma(link = "log"))
priors_default %>%
    gt()

prior	class	coef	group	resp	dpar	nlpar	lb	ub	source
	b								default
	b	cycles_z							default
	b	mass_z							default
	b	salinity_group9							default
	b	temp_mean_z							default
student_t(3, -2.8, 2.5)	Intercept								default
	b				shape				default
	b	salinity_group9			shape				default
student_t(3, 0, 2.5)	Intercept				shape				default

10.2.4 Run model

Here we run the model, I have hashed this out because I have saved the model for quick reloading.

⏭️ Skip this step if you have downloaded the saved model ‘smr_mod_gamma.rds’

Here we reload the model

setwd(output_mods_wd)

smr_mod_gamma <- readRDS(file = "smr_mod_gamma.rds")

10.2.5 Model diagnostics

Checking model convergence

color_scheme_set("red")
plot(smr_mod_gamma, ask = F)

Using leave one out (loo) measure of fit, the model appears to preform well, one Pareto k estimates falls outside the good range (0.7, 1]

loo(smr_mod_gamma)

## 
## Computed from 14000 by 58 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo    134.4  9.4
## p_loo         9.5  3.5
## looic      -268.9 18.9
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 0.9]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     56    96.6%   297     
##    (0.7, 1]   (bad)       2     3.4%   <NA>    
##    (1, Inf)   (very bad)  0     0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.

Model predictions generally align with the observed data, but there is a lot of uncertainty.

color_scheme_set("red")
p <- pp_check(smr_mod_gamma, type = "dens_overlay")
p

10.2.6 📈 Results

We did not see a meaningful difference between the routine metabolic rate for fish from the two salinity treatments.

10.2.6.1 Table S2

Table S2: Fixed effect Estimates (β) and 95% Credible Intervals (95% CI) from a Bayesian Generalised Linear Mixed Model (GLMM) with a Gamma distribution and a log link

model_est <- fixef(smr_mod_gamma, probs = c(0.025, 0.975)) %>%
    as.data.frame() %>%
    tibble::rownames_to_column(var = "Predictor") %>%
    dplyr::mutate(β = round(Estimate, 3), Q2.5 = round(Q2.5, 3), Q97.5 = round(Q97.5,
        3), `95% CI` = paste0("[", Q2.5, ", ", Q97.5, "]"))

model_est %>%
    dplyr::select(Predictor, "β", "95% CI") %>%
    gt()

Predictor	β	95% CI
Intercept	-2.751	[-2.889, -2.615]
shape_Intercept	2.018	[1.435, 2.545]
temp_mean_z	-0.364	[-1.548, 0.802]
cycles_z	-0.003	[-0.027, 0.019]
mass_z	1.379	[0.947, 1.818]
salinity_group9	-0.096	[-0.302, 0.11]
shape_salinity_group9	-0.081	[-0.913, 0.756]

Looking at the marginal mean difference between salinity groups

em_results <- emmeans(smr_mod_gamma, ~salinity_group)
contrast_results <- contrast(em_results, method = "pairwise")
em_results_df <- em_results %>%
    tidy() %>%
    mutate(across(where(is.numeric), ~exp(.)))
contrast_results_df <- contrast_results %>%
    tidy() %>%
    mutate(across(where(is.numeric), ~exp(.)))


em_results_df %>%
    gt()

salinity_group	estimate	lower.HPD	upper.HPD
0	0.06382226	0.05568229	0.07322901
9	0.05801414	0.04982011	0.06710878

Pulling the emmeans draws for our plot

emmeans_draws <- smr_mod_gamma %>%
    emmeans(~salinity_group) %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value), salinity_group = as.character(salinity_group))

emmeans_contrast_draws <- smr_mod_gamma %>%
    emmeans(~salinity_group) %>%
    contrast(method = "pairwise") %>%
    gather_emmeans_draws() %>%
    dplyr::mutate(.value = exp(.value))

10.2.6.2 Figure S9

This plot is in the main text of the manuscript as Figure 1b

mean_smr_salinity <- smr %>%
  dplyr::group_by(salinity_group) %>%
  dplyr::reframe(mean_SMR = mean(SMR, na.rm = TRUE))


fig_1b <- ggplot() +
  geom_violin(data = smr,
              aes(x = salinity_group, y = SMR, fill = salinity_group),
              color = NA, alpha = 0.2) +
  geom_jitter(data = smr,
              aes(x = salinity_group, y = SMR, fill = salinity_group),
              shape = 21, width = 0.3, size = 1, color = "black", alpha = 0.1) +
    geom_point(data = mean_smr_salinity,
             aes(x = salinity_group, y = mean_SMR, fill = salinity_group),
             size = 4, alpha = 1, stroke = 2, color = "black", shape = 21,
              position = position_nudge(x = 0.05)) +
  stat_pointinterval(data = emmeans_draws,
                     aes(x = salinity_group, y = .value),
                     color = "black", fill = "grey", point_interval = "mean_qi", .width = 0.95, shape = 21,  stroke = 2, point_size = 4, alpha = 1,
                     position = position_nudge(x = -0.05)) +
  scale_fill_manual(values = c("#4B5320", "#000080")) +  # Custom fill colours
  scale_colour_manual(values = c("#4B5320", "#000080")) +
  theme_clean() +
  theme(legend.position = "none") +
  labs(
    subtitle = "",
    x = "Salinity group (ppt)",
    y = "Standard metabolic rate (SMR; mg O2 g/h)"
  )

fig_1b

Figure 1b: The standard metabolic rate estimate (mg^-1 O₂ h^-1) plotted by salinity treatment. The small transparent points are the observed values, the shaded area behind the points is the a kernel density of the observed data, the large coloured point (to the right) is the observed mean, the large grey point with error bars (to the left) is the model estimated marginal mean (eemean) 95% Credible Intervals (95% CI).

10.3 Incremental regression analyses

Here we are following the methods Urbina et al. (2012)^[1] with an incremental regression analyses, in order to determine the best fit for the MO₂ vs O₂ data

This analysis approach evaluates the relative ‘fit’ of each polynomial order equation starting at zero and increasing to the third order, permitting a mathematical assessment of whether fish were oxyconforming or oxyregulatoring. If the data is best fitted/predicted by a single linear relationship (1^st-order polynomial) with a positive slope, this would suggests the fish were oxyconforming. Alternately, if the relationship is best modelled by a flat regression (0^th-order polynomial), or a higher order polynomial (2^nd or 3^rd-order polynomial) the fish is likely oxyregulatoring.

output_mods_bayes_wd <- paste0(output_mods_wd, "./bayes-regs")
ifelse(!dir.exists(output_mods_bayes_wd), dir.create(output_mods_bayes_wd), "Folder already exists")

## [1] "Folder already exists"

10.3.1 Building Bayesian regressions

Here we are using a Bayesian approach to model fitting with brm. These models take a long time to run, so I have saved them and re-loaded them to save time. I have also saved the summary data produced from the models, to save time, you can simply skip the hashed code and input the resulting summary data.

We will run our custom function, bayes_incremental_regression_by_id.

⏭️ Skip this step if you have already run this once, or have downloaded the saved models or saved datafiles from GitHub (that’s why its hashed out). This code takes a while to run

I have fund that I get an error, ‘Caused by error in socketConnection():’. I think the system may be hitting a limit on the number of simultaneous socket connections. If you try and run this code and also get this issue, try reducing the number of parallel workers: plan(multisession, workers = 2).

# ids <- slope_tidy %>% dplyr::distinct(id) %>% pull(id) plan(multisession)
# future_map( ids, bayes_incremental_regression_by_id, id_name = 'id', data =
# slope_tidy, response = 'MO2_g', predictor = 'DO', seed_number = 143019,
# save_models = TRUE, mod_output_wd = output_mods_bayes_wd ) plan(sequential)

Load all models and store in a list, will use a lot of memory. You can also skip this step and load the resulting data frames below. I am using the custom function load_rds, so we can compare them and generate predictions.

⏭️ Skip this step if you have downloaded the saved data pulled from these models have also hased this out, bayes_reg_mods_fit.csv’ and ‘bayes_reg_mods_predictions.csv’. These are in the mod-data directory.

# bayes_reg_mods <- load_rds(model_dw = output_mods_bayes_wd)

10.3.2 Model fits

⏭️ Skip this step if you have downloaded 💿 bayes_reg_mods_fit.csv.

It gets model fit parameters loo and r2 using the custom function, incremental_regression_bayes_fits. This code takes a while to run

# setwd(mod_data_wd) bayes_reg_mods_fit <-
# incremental_regression_bayes_fits(models = bayes_reg_mods)
# write.csv(bayes_reg_mods_fit, 'bayes_reg_mods_fit.csv', row.names = FALSE)

Reading in this model fit data frame, in the case you did not load in all the models.

setwd(mod_data_wd)
bayes_reg_mods_fit <- read.csv("bayes_reg_mods_fit.csv")

Selecting the best fitting model

elpd_loo, or the expected log pointwise predictive density for leave-one-out cross-validation, is a metric used in Bayesian model evaluation to assess the predictive accuracy of a model. The elpd_loo is an approximation of how well the model is expected to predict new data, based on leave-one-out cross-validation. Higher elpd_loo values indicate better predictive performance.

best_fit_bayes_reg <- bayes_reg_mods_fit %>%
    dplyr::group_by(id) %>%
    dplyr::mutate(elpd_loo_rank = rank(-elpd_loo)) %>%
    dplyr::select(id, model_type, elpd_loo, r2, elpd_loo_rank, r2_q2.5, r2_q97.5,
        estimate_DO, conf.low_DO, conf.high_DO) %>%
    ungroup()

10.3.3 Model predictions

⏭️ Skip this step if you have downloaded 💿 bayes_reg_mods_predictions.csv.

It pulls our model predictions using a custom function bayes_mod_predictions.

# setwd(mod_data_wd) bayes_reg_mods_predictions <- bayes_mod_predictions(models
# = bayes_reg_mods, original_data = slope_tidy)
# write.csv(bayes_reg_mods_predictions, 'bayes_reg_mods_predictions.csv',
# row.names = FALSE)

Reading in the predicted data

setwd(mod_data_wd)
bayes_reg_mods_predictions <- read.csv("bayes_reg_mods_predictions.csv")

We are going to combined this with our best fitting model df, so we know how they ranks for LOO.

bayes_reg_mods_predictions <- left_join(bayes_reg_mods_predictions, best_fit_bayes_reg,
    by = c("id", "model_type"))

10.3.4 📈 Results

10.3.4.1 Model selection summary

The best fitting models was most often a 2^nd-order polynomial (n = 25, 43%). This could suggest the presence of a critical oxygen threshold (Pcrit) where the relationship between O₂ and MO₂ changes. To confirm their is a Pcrit, we need to validated the shape of the polynomials and in should use a more specific model to test the Pcrit value. In any case, This type of model is indicative of oxyregulator. The same is also true for those with a 3^rd-order polynomial (n = 12, 21%).

The next most common is 1^st-order polynomial (n = 18, 31%). In the case of the 1^st-order polynomials, it suggest the presences of linear relationship between o2 and MO2, which is indicative of oxyconformer. However, to be true evidence of a oxyconformer this relationship should be positive (i.e. as O₂ falls MO₂ also falls). 13 of the 18 individuals best modelled with a linear function had positive estimates with credible intervals that did not overlap with zero (Table S3).

Lastly, 0^th-order was the least common (n = 3, 5%), it suggests that MO₂ does not show a statistically significant dependence on the O₂. In other words, the metabolic rate does not adjust based on oxygen availability, and there is no clear critical oxygen threshold (Pcrit) where the relationship changes. This is indicative of a oxyregulator.

best_mod <- best_fit_bayes_reg %>%
    dplyr::filter(elpd_loo_rank == 1)

total_fish <- nrow(best_mod)

mod_summary_table <- best_mod %>%
    dplyr::group_by(model_type) %>%
    dplyr::reframe(n = length(id), percent = round((n/total_fish) * 100, 2)) %>%
    dplyr::mutate(best_model_name = case_when(model_type == "lm_0" ~ "0th-order polynomial",
        model_type == "lm_1" ~ "1st-order polynomial", model_type == "lm_2" ~ "2nd-order polynomial",
        model_type == "lm_3" ~ "3rd-order polynomial", TRUE ~ "ERROR"))


table_bwm <- mod_summary_table %>%
    dplyr::select(best_model_name, everything(), -model_type) %>%
    gt() %>%
    cols_align(align = "center", columns = everything())

table_bwm

best_model_name	n	percent
0th-order polynomial	3	5.17
1st-order polynomial	18	31.03
2nd-order polynomial	25	43.10
3rd-order polynomial	12	20.69

Summary of fish best model with a linear function.

10.3.5 Table S3

Table S3: Ten fish that were best modelled with a linear function, showing r², and estimate (). Only two fish had positive estimates with credible intervals that did not overlap with zero, which are highlighted as conforming in the table. Thus in total, 13 of 58 fish showed sufficient evidence to be conforming.

table_lm_1 <- best_mod %>%
    dplyr::filter(model_type == "lm_1") %>%
    dplyr::mutate(r_sq_ci = paste0(round(r2, 3), " (", round(r2_q2.5, 3), "–",
        round(r2_q97.5, 3), ")"), est_ci = paste0(round(estimate_DO, 6), " (", round(conf.low_DO,
        6), "–", round(conf.high_DO, 6), ")"), conformer = if_else(conf.low_DO >
        0, "Conforming", "Not conforming")) %>%
    dplyr::select(id, r_sq_ci, est_ci, conformer) %>%
    dplyr::arrange(conformer) %>%
    dplyr::ungroup()


table_lm_1 %>%
    gt() %>%
    cols_align(align = "center", columns = everything()) %>%
    cols_label(id = "Fish ID", r_sq_ci = "r2 (CI)", est_ci = "Estimate (CI)", conformer = "Evidence of oxyconforming")

Fish ID	r2 (CI)	Estimate (CI)	Evidence of oxyconforming
a_0_24nov_1	0.187 (0.022–0.373)	0.000998 (0.000331–0.001701)	Conforming
a_0_24nov_3	0.254 (0.058–0.428)	0.001358 (0.000601–0.002094)	Conforming
a_0_25nov_1	0.537 (0.32–0.658)	0.000881 (0.000581–0.001173)	Conforming
a_0_25nov_4	0.253 (0.049–0.441)	0.000847 (0.000346–0.001356)	Conforming
a_9_22nov_1	0.198 (0.009–0.399)	0.00062 (0.000124–0.001096)	Conforming
b_0_25nov_2	0.451 (0.223–0.595)	0.000856 (0.000523–0.001177)	Conforming
b_0_27nov_3	0.254 (0.024–0.462)	0.000362 (0.000104–0.000614)	Conforming
b_9_22nov_1	0.421 (0.23–0.556)	0.000855 (0.00056–0.001146)	Conforming
c_9_24nov_4	0.294 (0.102–0.459)	0.000422 (0.000227–0.000623)	Conforming
c_9_25nov_4	0.535 (0.335–0.655)	0.000455 (0.00031–0.000601)	Conforming
c_9_27nov_2	0.186 (0.003–0.411)	0.00062 (7e-06–0.001206)	Conforming
d_0_22nov_3	0.141 (0.002–0.33)	0.000356 (2.6e-05–0.000677)	Conforming
d_9_25nov_3	0.294 (0.053–0.489)	0.001461 (0.000569–0.002336)	Conforming
a_9_21nov_1	0.131 (0.001–0.338)	0.000314 (-4.3e-05–0.000668)	Not conforming
b_0_24nov_4	0.101 (0.001–0.271)	0.000601 (-4e-05–0.001244)	Not conforming
b_0_27nov_2	0.136 (0.002–0.348)	0.000235 (-1.5e-05–0.000492)	Not conforming
c_0_22nov_3	0.109 (0–0.323)	0.000762 (-0.000341–0.001872)	Not conforming
d_9_25nov_2	0.124 (0.001–0.33)	0.000386 (-5.9e-05–0.000828)	Not conforming

10.3.5.1 Ploting all models

Data set with all slopes and which model was best

best_fit <- left_join(slope_tidy, best_mod, by = "id")

Now we are plotting each of the regressions. First making a directory to save the figures

Plot all regression, and highlighting the model that has the best fit, based on AIC values

# Create a list to store the plots
plots <- list()
model_preds_list <- list()

for (id_i in ids) {
  
  # Filter data for the current ID
  df_i <- bayes_reg_mods_predictions %>%
    dplyr::filter(id == id_i) %>% 
    dplyr::mutate(line_size = if_else(elpd_loo_rank == 1, 2, 1),
           alpha_value = if_else(elpd_loo_rank == 1, 1, 0.4))
  
  x_min <- df_i %>%
    dplyr::reframe(min = min(DO), na.rm = TRUE) %>% 
    dplyr::pull(min)
  
  y_max <- df_i %>%
    dplyr::reframe(max = max(MO2_g), na.rm = TRUE) %>% 
    dplyr::pull(max)
  
  best_weighted_model_i <- best_fit_bayes_reg %>% 
    dplyr::filter(id == id_i & elpd_loo_rank == 1)
  
  poly_i_name <- best_weighted_model_i %>%
    dplyr::mutate(name = case_when(
      model_type == "lm_0" ~ "0th-order",
      model_type == "lm_1" ~ "1st-order",
      model_type == "lm_2" ~ "2nd-order",
      model_type == "lm_3" ~ "3rd-order",
      TRUE ~ "ERROR"
    )) %>% 
    dplyr::pull(name)
  
  r_i <- best_weighted_model_i %>% 
    dplyr::mutate(r_sq_ci = paste0(round(r2, 3), " (", 
                                    round(r2_q2.5, 3), "–", 
                                    round(r2_q97.5, 3), ")")) %>% 
    dplyr::pull(r_sq_ci)

  # Create the plot
  plot_i <- ggplot() +
    geom_ribbon(data = df_i,
                aes(x = DO, y = predicted, ymin = pred_lower, ymax = pred_upper, fill = model_type),
                alpha = 0.1) +
    geom_line(data = df_i, 
              aes(x = DO, y = predicted, colour = model_type, size = line_size, alpha = alpha_value)) +
    geom_point(data = df_i %>% dplyr::filter(elpd_loo_rank == 1), aes(x = DO, y = MO2_g), alpha = 0.6, colour = "black", size = 2) +
    scale_colour_manual(values = c("red", "blue", "green", "purple"), 
                        labels = c("0th Order", "1st Order", "2nd Order", "3rd Order")) +
    scale_size_identity() +  # Use the size values directly
    scale_alpha_identity(guide = "none") +  # Remove the alpha legend 
    annotate("text", x = x_min, 
             y = y_max, 
             label = paste0("Best fit: ",poly_i_name, "\n", "r2 = ", r_i), 
             hjust = 0, vjust = 1, size = 4) +
    labs(
      title = paste(id_i),
      x = "Dissolved oxygen percentage (DO)",
      y = "MO2 (o2 mg/g/h)",
      colour = "Model") +
    theme_classic() +
    theme(legend.position = "none")
  
  # Store the plot
  plots[[id_i]] <- plot_i
  
  print(plot_i)
}

pdf(file = paste0(incremental_reg_bayes_wd, "./combined_reg_plots.pdf"), width = 8, height = 6)

for (id_i in ids) {
  print(plots[[id_i]])  # Print each plot to the PDF
}

dev.off()  # Close the PDF device

## png 
##   2

Getting a plot for each best fit regression, and overlaying a global model for that model type.

## [1] "Folder already exists"

Here we are grouping fish by best fitting model and getting an average trend. I have hashed the code so you don’t need to re-run the models

⏭️ Skip this step if you have downloaded the global models ‘.rds’ in the ‘bayes-regs-global’ folder

# setwd(output_mods_bayes_global_wd) seed_number = 143019 # Set up parallel
# backend plan(multisession, workers = 4) # Adjust workers based on system
# resources # Define model formulas and data model_formulas <- list( lm_0 =
# bf(MO2_g ~ 1, family = gaussian()), lm_1 = bf(MO2_g ~ DO, family =
# gaussian()), lm_2 = bf(MO2_g ~ poly(DO, 2), family = gaussian()), lm_3 =
# bf(MO2_g ~ poly(DO, 3), family = gaussian()) ) model_data <- list( lm_0 =
# best_fit %>% filter(model_type == 'lm_0'), lm_1 = best_fit %>%
# filter(model_type == 'lm_1'), lm_2 = best_fit %>% filter(model_type ==
# 'lm_2'), lm_3 = best_fit %>% filter(model_type == 'lm_3') ) # Run models in
# parallel with future_map() models <- future_map(names(model_formulas), ~{
# brm( formula = model_formulas[[.x]], data = model_data[[.x]], cores = 4, seed
# = seed_number, save_pars = save_pars(all = TRUE), sample_prior = FALSE,
# silent = TRUE, file = paste0(.x, '_global') ) }) # Assign model names to
# results names(models) <- names(model_formulas) # Stop parallel plan after
# execution plan(sequential)

Loading these models

bayes_reg_mods_global <- load_rds(model_dw = output_mods_bayes_global_wd)

Pulling predictions

global_models_preds <- list()

for (mode_name in names(bayes_reg_mods_global)) {

    mod_i <- bayes_reg_mods_global[[mode_name]]
    mode_type = str_remove(mode_name, "_global")

    model_predictions_i <- as.data.frame(fitted(mod_i, summary = TRUE)) %>%
        dplyr::mutate(model = mode_name, model_type = mode_type) %>%
        dplyr::rename(pred_lower = Q2.5, pred_upper = Q97.5, predicted = Estimate,
            pred_error = Est.Error) %>%
        dplyr::select(model, everything())

    original_data_i <- best_fit %>%
        filter(model_type == mode_type) %>%
        dplyr::select(-model_type)

    model_predictions_original_i <- cbind(model_predictions_i, original_data_i)

    global_models_preds[[mode_name]] <- model_predictions_original_i
}

global_models_pred_df <- bind_rows(global_models_preds)

10.3.5.2 Figure S10

Best fit regressions with global models based on that polynomial order.

Figure S10 is also called Figure 2 in the manuscript.

global_models_pred_df <- global_models_pred_df %>%
  dplyr::mutate(model_name = case_when(
      model_type == "lm_0" ~ "0th-order polynomial",
      model_type == "lm_1" ~ "1st-order polynomial",
      model_type == "lm_2" ~ "2nd-order polynomial",
      model_type == "lm_3" ~ "3rd-order polynomial",
      TRUE ~ "ERROR"
    ))

annotation_data <- mod_summary_table %>%
  dplyr::select(model_type, best_model_name, n)

bayes_reg_mods_predictions_best <- bayes_reg_mods_predictions %>% 
  dplyr::filter(elpd_loo_rank == 1) 

fig_1 <- ggplot() +
  geom_line(data = bayes_reg_mods_predictions_best,
            aes(x = DO, y = predicted, color = id), size = 1) +
  # geom_smooth(data = best_fit,
  #           aes(x = DO, y = MO2_g, color = id), size = 1, alpha = 1, se = FALSE) +
  geom_ribbon(data = global_models_pred_df,
              aes(x = DO, ymin = pred_lower
, ymax = pred_upper, group = model),
              fill = "#FC6C85", alpha = 0.2) +  # Shaded confidence intervals
  geom_line(data = global_models_pred_df,
            aes(x = DO, y = predicted), size = 1.5, color = "#FF007F") +
  facet_wrap(~model_type) +
  scale_color_grey(start = 0.1, end = 0.9) +
  labs(
      title = paste("Best fit regressions grouped by polynomial order"),
      x = "Dissolved oxygen percentage (DO)",
      y = "MO2 (O2 mg/g/h)") +
  theme_classic() +
  theme(legend.position = "none") +
  geom_text(data = annotation_data,
            aes(x = -Inf, y = Inf, label = paste0("italic(n) == ", n)),
            hjust = -0.1, vjust = 1.2, inherit.aes = FALSE, parse = TRUE)
fig_1

Figure 2: Best fit models from incremental regression analyses (0^th to 3 ^rd polynomial regressions) grouped by polynomial order, the pink regression represents the average for all fish fit with that polynomial order.

10.4 Pcrit model

For those fish that were best modelled with a 2^nd or 3^rd-order polynomial (n = 38) we will check to see if a Pcrit is present. We are filtering the data for only those fish.

higher_order <- c("lm_2", "lm_3")

check_pcrit_ids <- best_mod %>%
    dplyr::filter(model_type %in% higher_order) %>%
    dplyr::distinct(id) %>%
    pull(id)

check_pcrit_df <- slope_tidy %>%
    dplyr::filter(id %in% check_pcrit_ids)

We will calculate Pcrit using Chabot method and function calcO2crit.

This function uses the fifth percentile of the MO₂ values observed at dissolved oxygen levels ≥ 80% air saturation as the criterion to assess low MO₂ values. The algorithm then identifies all the MO₂ measurements greater than this minimally acceptable MO₂ value. Within this sub-set, it identifies the ̇ MO2 measurement made at the lowest DO and thereafter considers this DO as candidate for breakpoint (named pivotDO in the script). A regression is then calculated using observations at DO levels < pivotDO, and a first estimate of O2crit is calculated as the intersection of this regression line with the horizontal line representing SMR. The script then goes through validation steps to ensure that the slope of the regression is not so low that the line, projected to normoxic DO levels, passes below any MO2 values observed in normoxia. It also ensures that the intercept is not greater than zero. Corrective measures are taken if such problems are encountered.

lowestMO2 default is the quantile(MO2[DO >= 80], p=0.05). It is used to segment the data and locate the pivotDO.

ids <- check_pcrit_df %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

pcrit_model_df_list <- list()
pcrit_models <- list()

for (id_i in ids) {

    df_i <- check_pcrit_df %>%
        dplyr::filter(id == id_i)

    o2crit <- calcO2crit(Data = df_i, SMR = df_i$SMR[1], lowestMO2 = NA, gapLimit = 4,
        max.nb.MO2.for.reg = 7)

    vaule <- o2crit$o2crit
    lowestMO2 = quantile(df_i$MO2[df_i$DO >= 80], p = 0.05)
    SMR <- o2crit$SMR
    nb_mo2_conforming <- o2crit$Nb_MO2_conforming
    r2 <- o2crit$r2
    method <- o2crit$Method
    p <- o2crit$P[1]

    pcrit_model_df <- tibble(id = id_i, pcrit_vaule = vaule, pcrit_smr = SMR, pcrit_lowestMO2 = lowestMO2,
        pcrit_nb_mo2_conforming = nb_mo2_conforming, pcrit_r2 = r2, pcrit_method = method,
        pcrit_p = p)

    pcrit_model_df_list[[id_i]] <- pcrit_model_df

    pcrit_models[[id_i]] <- o2crit

}

pcrit_model_df <- bind_rows(pcrit_model_df_list)

10.4.1 Ploting Pcrit

Here’s the plots for the Pcrit estimates. This inculdes all fish with best fit regression that are 2 or 3 order (n = 38).

## png 
##   2

Printing in htlm document

ids <- check_pcrit_df %>%
    dplyr::distinct(id) %>%
    dplyr::pull()

for (id_i in ids) {

    comment <- check_pcrit_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::slice(1) %>%
        dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
        pull(comment)

    r2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
        dplyr::pull(pcrit_r2)

    conforming <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
        dplyr::pull(pcrit_nb_mo2_conforming)

    P <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
        dplyr::pull(pcrit_p)

    SMR <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
        dplyr::pull(pcrit_smr)

    lowestMO2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
        dplyr::pull(pcrit_lowestMO2)

    # Generate and render the plot
    plotO2crit(o2critobj = pcrit_models[[id_i]])

    # Add a title
    mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
        font = 2, cex = 1.2)

    mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
        SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
        font = 0.5, cex = 0.8)
}

10.4.2 🚧 Visual inspection 🚧

🚧 Need to re-do visual inspections with new data 🚧

We conducted visual inspection of all 38 fish to detrmine if Pcrits were present. We have laoded this as a data frame 💿 pcrit_check

For this visual inspection, we grouped fish into “no”, “unclear” or “yes”, based on the certainty at which we observed a Pcrit. We concluded that 12 fish had a Pcrit responses, while 6 showed some evidence of a Pcrit but not as clearly (i.e. unclear), and the remaining 20 did not show any evidence of a Pcrit response.

pcrit_check %>%
    dplyr::group_by(pcrit) %>%
    dplyr::reframe(n = length(id)) %>%
    gt()

pcrit	n
no	20
unclear	6
yes	12

Let’s make a list of those that we did visually identify as a Pcrit

pcrit_visual_yes <- pcrit_check %>%
    dplyr::filter(pcrit == "yes") %>%
    dplyr::pull(id)

10.4.3 Pcrit numerical rules

We also used a numerical rule to asses if a Pcrit was observed.

We filtered for only cases were at the lowest O₂ value three consecutive MO₂ measures full below our SMR and fifth percentile of the MO2 values observed at dissolved O2 levels > 80%. In the model output these are called nb_mo2_conforming points.

pcrit_numerical_yes <- pcrit_model_df %>%
    dplyr::filter(pcrit_nb_mo2_conforming > 2) %>%
    pull(id)

paste0("Based on this rule there are ", length(pcrit_numerical_yes), " fish with possible Pcrits.")

## [1] "Based on this rule there are 21 fish with possible Pcrits."

10.4.4 🚧 📈 Pcrits numbers 🚧

🚧 Need to re-do visual inspections with new data 🚧

We will now check to see which fish have Pcrits based on the numerical and visual inspections.

shared_elements <- intersect(pcrit_visual_yes, pcrit_numerical_yes)

unique_to_visual <- setdiff(pcrit_visual_yes, pcrit_numerical_yes)

unique_to_numerical <- setdiff(pcrit_numerical_yes, pcrit_visual_yes)

list(Shared = shared_elements, Unique_to_visual = unique_to_visual, Unique_to_numerical = unique_to_numerical)

## $Shared
##  [1] "a_0_26nov_4" "a_9_21nov_3" "b_0_24nov_1" "b_0_24nov_2" "b_0_25nov_1"
##  [6] "b_0_25nov_3" "b_0_26nov_1" "b_9_21nov_1" "b_9_21nov_2" "b_9_21nov_3"
## [11] "d_0_21nov_3"
## 
## $Unique_to_visual
## [1] "d_0_21nov_2"
## 
## $Unique_to_numerical
##  [1] "a_0_24nov_4" "b_0_25nov_4" "b_0_26nov_3" "b_9_21nov_4" "b_9_22nov_3"
##  [6] "c_0_21nov_2" "c_9_24nov_2" "c_9_27nov_4" "d_0_22nov_2" "d_9_24nov_3"

There are 10 fish that both numerically and visually were determined to have a Pcrit. Two fish were detriment visually to have a Pcrit, while they did not meet the numerical criteria. Three fish meet the numerical criteria, but did not pass the visual inspection.

These are the 10 fish that meet numerical criteria and passed visual inspection.

for (id_i in shared_elements) {

    comment <- check_pcrit_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::slice(1) %>%
        dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
        pull(comment)

    r2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
        dplyr::pull(pcrit_r2)

    conforming <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
        dplyr::pull(pcrit_nb_mo2_conforming)

    P <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
        dplyr::pull(pcrit_p)

    SMR <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
        dplyr::pull(pcrit_smr)

    lowestMO2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
        dplyr::pull(pcrit_lowestMO2)

    # Generate and render the plot
    plotO2crit(o2critobj = pcrit_models[[id_i]])

    # Add a title
    mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
        font = 2, cex = 1.2)

    mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
        SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
        font = 1, cex = 0.8)
}

These are the 2 fish that did not meet numerical criteria but passed visual inspection

for (id_i in unique_to_visual) {

    comment <- check_pcrit_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::slice(1) %>%
        dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
        pull(comment)

    r2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
        dplyr::pull(pcrit_r2)

    conforming <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
        dplyr::pull(pcrit_nb_mo2_conforming)

    P <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
        dplyr::pull(pcrit_p)

    SMR <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
        dplyr::pull(pcrit_smr)

    lowestMO2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
        dplyr::pull(pcrit_lowestMO2)

    # Generate and render the plot
    plotO2crit(o2critobj = pcrit_models[[id_i]])

    # Add a title
    mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
        font = 2, cex = 1.2)

    mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
        SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
        font = 1, cex = 0.8)
}

These are the 3 fish that did meet numerical criteria but did not pass visual inspection

for (id_i in unique_to_numerical) {

    comment <- check_pcrit_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::slice(1) %>%
        dplyr::mutate(comment = if_else(is.na(comments), "", paste0("#", comments))) %>%
        pull(comment)

    r2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_r2 = round(pcrit_r2, 3)) %>%
        dplyr::pull(pcrit_r2)

    conforming <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_nb_mo2_conforming = round(pcrit_nb_mo2_conforming, 3)) %>%
        dplyr::pull(pcrit_nb_mo2_conforming)

    P <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_p = round(pcrit_p, 3)) %>%
        dplyr::pull(pcrit_p)

    SMR <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_smr = round(pcrit_smr, 3)) %>%
        dplyr::pull(pcrit_smr)

    lowestMO2 <- pcrit_model_df %>%
        dplyr::filter(id == id_i) %>%
        dplyr::mutate(pcrit_lowestMO2 = round(pcrit_lowestMO2, 3)) %>%
        dplyr::pull(pcrit_lowestMO2)

    # Generate and render the plot
    plotO2crit(o2critobj = pcrit_models[[id_i]])

    # Add a title
    mtext(text = paste0(id_i, " ", comment), side = 3, line = 2, adj = 0, col = "blue",
        font = 2, cex = 1.2)

    mtext(text = paste0("R2 = ", r2, "; p = ", P, "; CP < SMR = ", conforming, "; SMR = ",
        SMR, "; lowestMO2 = ", lowestMO2), side = 3, line = 1, adj = 0, col = "blue",
        font = 1, cex = 0.8)
}

Based on the 10 fish that passed both the numerical criteria and visually checks the Pcrit response.

n_pcrit <- length(shared_elements)

have_pcirt <- pcrit_model_df %>%
    dplyr::filter(id %in% shared_elements)

mean_pcrit <- have_pcirt %>%
    dplyr::reframe(mean = mean(pcrit_vaule)) %>%
    pull(mean) %>%
    round(., 2)

min_pcrit <- have_pcirt %>%
    dplyr::reframe(min = min(pcrit_vaule)) %>%
    pull(min) %>%
    round(., 2)

max_pcrit <- have_pcirt %>%
    dplyr::reframe(max = max(pcrit_vaule)) %>%
    pull(max) %>%
    round(., 2)

print(paste0("There are ", n_pcrit, " fish with identified Pcrits and the mean Pcrit is ",
    mean_pcrit, " (range: ", min_pcrit, "-", max_pcrit, ")", " however, the true range for Pcrit is likey much larger then reported here, as many fish did not exhibit a Pcrit like response with the trial"))

## [1] "There are 11 fish with identified Pcrits and the mean Pcrit is 32.78 (range: 21.3-57.4) however, the true range for Pcrit is likey much larger then reported here, as many fish did not exhibit a Pcrit like response with the trial"

The mean lowest dissolved oxygen (DO) value for trails where a Pcrit was not detected was 13.11% ± 3.30% (mean ± sd), the lowest recoded DO value for which a Pcrit was not detected was 8.82%. Thus the true Pcrit range of Galaxias maculatus is 39.9% to < 8.83%

slope_tidy %>%
    dplyr::filter(!(id %in% shared_elements)) %>%
    dplyr::group_by(id) %>%
    dplyr::arrange(DO) %>%
    dplyr::slice(1) %>%
    dplyr::ungroup() %>%
    dplyr::reframe(mean_lowest_o2 = mean(DO), sd_lowest_o2 = sd(DO), min_o2 = min(DO))

## # A tibble: 1 × 3
##   mean_lowest_o2 sd_lowest_o2 min_o2
##            <dbl>        <dbl>  <dbl>
## 1           13.1         3.33   8.82

For the fish that did have Pcrits, proportionally more fish were from the zero salinity treatment (20% vs 14%), and the fish were close to the mean weight.

mean_weight_all <- slope_tidy %>%
    dplyr::filter(!(id %in% shared_elements)) %>%
    dplyr::group_by(id) %>%
    dplyr::slice(1) %>%
    dplyr::ungroup() %>%
    dplyr::reframe(mean_weight_all = mean(mass, na.rm = TRUE)) %>%
    dplyr::pull(mean_weight_all)

total_salinity <- slope_tidy %>%
    dplyr::group_by(id) %>%
    dplyr::slice(1) %>%
    dplyr::ungroup() %>%
    dplyr::group_by(salinity_group) %>%
    dplyr::reframe(n = length(id))

n_salinity_0 <- total_salinity %>%
    dplyr::filter(salinity_group == "0") %>%
    dplyr::pull(n)

n_salinity_9 <- total_salinity %>%
    dplyr::filter(salinity_group == "9") %>%
    dplyr::pull(n)

slope_tidy %>%
    dplyr::filter(id %in% shared_elements) %>%
    dplyr::group_by(id) %>%
    dplyr::slice(1) %>%
    dplyr::ungroup() %>%
    dplyr::reframe(n_pcrit_salinity_9 = sum(salinity_group == "9", na.rm = TRUE),
        n_pcrit_salinity_0 = sum(salinity_group == "0", na.rm = TRUE), prop_salinity_9 = n_pcrit_salinity_9/n_salinity_9,
        prop_salinity_0 = n_pcrit_salinity_0/n_salinity_0, rel_mean_weight = mean(mass,
            na.rm = TRUE)/mean_weight_all) %>%
    gt()

n_pcrit_salinity_9	n_pcrit_salinity_0	prop_salinity_9	prop_salinity_0	rel_mean_weight
4	7	0.1428571	0.2333333	0.903909

11 📑 References

^[1] Claireaux, G. and Chabot, D. (2016) Responses by fishes to environmental hypoxia: integration through Fry’s concept of aerobic metabolic scope. Journal of Fish Biology https://doi.org/10.1111/jfb.12833

^[2] Urbina MA, Glover CN, and Forster ME, (2012) A novel oxyconforming response in the freshwater fish Galaxias maculatus. Comparative Biochemistry and Physiology Part A: Molecular & Integrative Physiology. https://doi.org/10.1016/j.cbpa.2011.11.011

^[3] Pick JL, Nakagawa S, and Noble DWA (2018) Reproducible, flexible and high-throughput data extraction from primary literature: The metaDigitise r package. https://doi.org/10.1111/2041-210X.13118

12 💻 Session information

Here is a detailed list of the session information

sessioninfo::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.2.3 (2023-03-15 ucrt)
##  os       Windows 10 x64 (build 19044)
##  system   x86_64, mingw32
##  ui       RTerm
##  language (EN)
##  collate  English_Australia.utf8
##  ctype    English_Australia.utf8
##  tz       Australia/Sydney
##  date     2025-02-18
##  pandoc   3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package           * version    date (UTC) lib source
##  abind               1.4-5      2016-07-21 [1] CRAN (R 4.2.0)
##  arrayhelpers        1.1-0      2020-02-04 [1] CRAN (R 4.2.3)
##  backports           1.4.1      2021-12-13 [1] CRAN (R 4.2.0)
##  bayesplot         * 1.11.1     2024-02-15 [1] CRAN (R 4.2.3)
##  betareg           * 3.2-0      2024-07-07 [1] CRAN (R 4.2.3)
##  bitops              1.0-7      2021-04-24 [1] CRAN (R 4.2.0)
##  boot                1.3-28.1   2022-11-22 [2] CRAN (R 4.2.3)
##  bridgesampling      1.1-2      2021-04-16 [1] CRAN (R 4.2.3)
##  brms              * 2.21.0     2024-03-20 [1] CRAN (R 4.2.3)
##  Brobdingnag         1.2-9      2022-10-19 [1] CRAN (R 4.2.3)
##  broom             * 1.0.4      2023-03-11 [1] CRAN (R 4.2.3)
##  broom.helpers       1.15.0     2024-04-05 [1] CRAN (R 4.2.3)
##  broom.mixed       * 0.2.9.5    2024-04-01 [1] CRAN (R 4.2.3)
##  bslib               0.7.0      2024-03-29 [1] CRAN (R 4.2.3)
##  cachem              1.0.7      2023-02-24 [1] CRAN (R 4.2.3)
##  car               * 3.1-2      2023-03-30 [1] CRAN (R 4.2.3)
##  carData           * 3.0-5      2022-01-06 [1] CRAN (R 4.2.3)
##  caTools             1.18.2     2021-03-28 [1] CRAN (R 4.2.3)
##  cellranger          1.1.0      2016-07-27 [1] CRAN (R 4.2.3)
##  checkmate           2.3.1      2023-12-04 [1] CRAN (R 4.2.3)
##  cli                 3.6.1      2023-03-23 [1] CRAN (R 4.2.3)
##  cluster             2.1.4      2022-08-22 [2] CRAN (R 4.2.3)
##  coda                0.19-4.1   2024-01-31 [1] CRAN (R 4.2.3)
##  codetools           0.2-19     2023-02-01 [2] CRAN (R 4.2.3)
##  colorspace          2.1-0      2023-01-23 [1] CRAN (R 4.2.3)
##  corrplot          * 0.92       2021-11-18 [1] CRAN (R 4.2.3)
##  crosstalk           1.2.1      2023-11-23 [1] CRAN (R 4.2.3)
##  curl                5.2.1      2024-03-01 [1] CRAN (R 4.2.3)
##  dagitty           * 0.3-4      2023-12-07 [1] CRAN (R 4.2.3)
##  data.table        * 1.14.8     2023-02-17 [1] CRAN (R 4.2.3)
##  datawizard        * 0.12.0     2024-07-11 [1] CRAN (R 4.2.3)
##  DEoptimR            1.1-3      2023-10-07 [1] CRAN (R 4.2.3)
##  devtools          * 2.4.5      2022-10-11 [1] CRAN (R 4.2.3)
##  digest              0.6.31     2022-12-11 [1] CRAN (R 4.2.3)
##  distributional      0.4.0      2024-02-07 [1] CRAN (R 4.2.3)
##  doParallel          1.0.17     2022-02-07 [1] CRAN (R 4.2.3)
##  dplyr             * 1.1.1      2023-03-22 [1] CRAN (R 4.2.3)
##  DT                * 0.33       2024-04-04 [1] CRAN (R 4.2.3)
##  ellipsis            0.3.2      2021-04-29 [1] CRAN (R 4.2.3)
##  emmeans           * 1.10.3     2024-07-01 [1] CRAN (R 4.2.3)
##  estimability        1.5.1      2024-05-12 [1] CRAN (R 4.2.3)
##  evaluate            0.24.0     2024-06-10 [1] CRAN (R 4.2.3)
##  fansi               1.0.4      2023-01-22 [1] CRAN (R 4.2.3)
##  farver              2.1.1      2022-07-06 [1] CRAN (R 4.2.3)
##  fastmap             1.1.1      2023-02-24 [1] CRAN (R 4.2.3)
##  flexmix             2.3-19     2023-03-16 [1] CRAN (R 4.2.3)
##  fontawesome         0.5.2      2023-08-19 [1] CRAN (R 4.2.3)
##  forcats           * 1.0.0      2023-01-29 [1] CRAN (R 4.2.3)
##  foreach             1.5.2      2022-02-02 [1] CRAN (R 4.2.3)
##  formatR             1.14       2023-01-17 [1] CRAN (R 4.2.3)
##  Formula             1.2-5      2023-02-24 [1] CRAN (R 4.2.2)
##  fs                  1.6.4      2024-04-25 [1] CRAN (R 4.2.3)
##  furrr             * 0.3.1      2022-08-15 [1] CRAN (R 4.2.3)
##  future            * 1.33.2     2024-03-26 [1] CRAN (R 4.2.3)
##  generics            0.1.3      2022-07-05 [1] CRAN (R 4.2.3)
##  ggdag             * 0.2.12     2024-03-08 [1] CRAN (R 4.2.3)
##  ggdist              3.3.2      2024-03-05 [1] CRAN (R 4.2.3)
##  ggExtra           * 0.10.0     2022-03-23 [1] CRAN (R 4.2.3)
##  ggfortify         * 0.4.16     2023-03-20 [1] CRAN (R 4.2.3)
##  gghalves          * 0.1.4      2022-11-20 [1] CRAN (R 4.2.3)
##  ggplot2           * 3.5.1      2024-04-23 [1] CRAN (R 4.2.3)
##  ggridges          * 0.5.6      2024-01-23 [1] CRAN (R 4.2.3)
##  ggthemes          * 4.2.4      2021-01-20 [1] CRAN (R 4.2.3)
##  globals             0.16.3     2024-03-08 [1] CRAN (R 4.2.3)
##  glue                1.6.2      2022-02-24 [1] CRAN (R 4.2.3)
##  gridExtra         * 2.3        2017-09-09 [1] CRAN (R 4.2.3)
##  gsw                 1.2-0      2024-08-19 [1] CRAN (R 4.2.3)
##  gt                * 0.11.0     2024-07-09 [1] CRAN (R 4.2.3)
##  gtable              0.3.5      2024-04-22 [1] CRAN (R 4.2.3)
##  gtExtras          * 0.4.5      2022-11-28 [1] CRAN (R 4.2.3)
##  gtsummary         * 1.7.2      2023-07-15 [1] CRAN (R 4.2.3)
##  highr               0.11       2024-05-26 [1] CRAN (R 4.2.3)
##  hms               * 1.1.3      2023-03-21 [1] CRAN (R 4.2.3)
##  htmltools           0.5.8.1    2024-04-04 [1] CRAN (R 4.2.3)
##  htmlwidgets         1.6.2      2023-03-17 [1] CRAN (R 4.2.3)
##  httpuv              1.6.9      2023-02-14 [1] CRAN (R 4.2.3)
##  httr                1.4.7      2023-08-15 [1] CRAN (R 4.2.3)
##  igraph            * 1.6.0      2023-12-11 [1] CRAN (R 4.2.3)
##  inline              0.3.19     2021-05-31 [1] CRAN (R 4.2.3)
##  insight             0.20.1     2024-06-11 [1] CRAN (R 4.2.3)
##  iterators           1.0.14     2022-02-05 [1] CRAN (R 4.2.3)
##  janitor           * 2.2.0      2023-02-02 [1] CRAN (R 4.2.3)
##  jquerylib           0.1.4      2021-04-26 [1] CRAN (R 4.2.3)
##  jsonlite            1.8.4      2022-12-06 [1] CRAN (R 4.2.3)
##  knitr               1.42       2023-01-25 [1] CRAN (R 4.2.3)
##  labeling            0.4.3      2023-08-29 [1] CRAN (R 4.2.3)
##  later               1.3.0      2021-08-18 [1] CRAN (R 4.2.3)
##  lattice           * 0.20-45    2021-09-22 [2] CRAN (R 4.2.3)
##  lazyeval            0.2.2      2019-03-15 [1] CRAN (R 4.2.3)
##  lifecycle           1.0.4      2023-11-07 [1] CRAN (R 4.2.3)
##  listenv             0.9.1      2024-01-29 [1] CRAN (R 4.2.3)
##  lme4              * 1.1-35.5   2024-07-03 [1] CRAN (R 4.2.3)
##  lmerTest          * 3.1-3      2020-10-23 [1] CRAN (R 4.2.3)
##  lmtest              0.9-40     2022-03-21 [1] CRAN (R 4.2.3)
##  loo                 2.8.0      2024-07-03 [1] CRAN (R 4.2.3)
##  lubridate         * 1.9.2      2023-02-10 [1] CRAN (R 4.2.3)
##  magrittr            2.0.3      2022-03-30 [1] CRAN (R 4.2.3)
##  marelac             2.1.11     2023-09-25 [1] CRAN (R 4.2.3)
##  marginaleffects   * 0.21.0     2024-06-14 [1] CRAN (R 4.2.3)
##  MASS                7.3-58.2   2023-01-23 [2] CRAN (R 4.2.3)
##  Matrix            * 1.5-3      2022-11-11 [2] CRAN (R 4.2.3)
##  matrixStats         1.3.0      2024-04-11 [1] CRAN (R 4.2.3)
##  mclust            * 6.1.1      2024-04-29 [1] CRAN (R 4.2.3)
##  measurements        1.5.1      2023-04-28 [1] CRAN (R 4.2.3)
##  memoise             2.0.1      2021-11-26 [1] CRAN (R 4.2.3)
##  mgcv                1.9-1      2023-12-21 [1] CRAN (R 4.2.3)
##  mime                0.12       2021-09-28 [1] CRAN (R 4.2.0)
##  miniUI              0.1.1.1    2018-05-18 [1] CRAN (R 4.2.3)
##  minqa               1.2.7      2024-05-20 [1] CRAN (R 4.2.3)
##  modeltools          0.2-23     2020-03-05 [1] CRAN (R 4.2.0)
##  multcomp            1.4-25     2023-06-20 [1] CRAN (R 4.2.3)
##  munsell             0.5.1      2024-04-01 [1] CRAN (R 4.2.3)
##  mvtnorm             1.2-5      2024-05-21 [1] CRAN (R 4.2.3)
##  nlme                3.1-162    2023-01-31 [2] CRAN (R 4.2.3)
##  nloptr              2.1.1      2024-06-25 [1] CRAN (R 4.2.3)
##  nnet                7.3-18     2022-09-28 [2] CRAN (R 4.2.3)
##  numDeriv            2016.8-1.1 2019-06-06 [1] CRAN (R 4.2.0)
##  oce                 1.8-3      2024-08-17 [1] CRAN (R 4.2.3)
##  opdisDownsampling   1.0.1      2024-04-15 [1] CRAN (R 4.2.3)
##  paletteer           1.6.0      2024-01-21 [1] CRAN (R 4.2.3)
##  parallelly          1.37.1     2024-02-29 [1] CRAN (R 4.2.3)
##  pbmcapply           1.5.1      2022-04-28 [1] CRAN (R 4.2.0)
##  performance       * 0.12.0     2024-06-08 [1] CRAN (R 4.2.3)
##  permute           * 0.9-7      2022-01-27 [1] CRAN (R 4.2.3)
##  pillar              1.9.0      2023-03-22 [1] CRAN (R 4.2.3)
##  pkgbuild            1.4.4      2024-03-17 [1] CRAN (R 4.2.3)
##  pkgconfig           2.0.3      2019-09-22 [1] CRAN (R 4.2.3)
##  pkgload             1.3.2      2022-11-16 [1] CRAN (R 4.2.3)
##  plotly            * 4.10.1     2022-11-07 [1] CRAN (R 4.2.3)
##  plyr                1.8.8      2022-11-11 [1] CRAN (R 4.2.3)
##  posterior           1.6.0      2024-07-03 [1] CRAN (R 4.2.3)
##  pracma              2.4.4      2023-11-10 [1] CRAN (R 4.2.3)
##  profvis             0.3.8      2023-05-02 [1] CRAN (R 4.2.3)
##  promises            1.2.0.1    2021-02-11 [1] CRAN (R 4.2.3)
##  purrr             * 1.0.1      2023-01-10 [1] CRAN (R 4.2.3)
##  qqconf              1.3.2      2023-04-14 [1] CRAN (R 4.2.3)
##  qqplotr           * 0.0.6      2023-01-25 [1] CRAN (R 4.2.3)
##  QuickJSR            1.3.0      2024-07-08 [1] CRAN (R 4.2.3)
##  R6                  2.5.1      2021-08-19 [1] CRAN (R 4.2.3)
##  RColorBrewer      * 1.1-3      2022-04-03 [1] CRAN (R 4.2.0)
##  Rcpp              * 1.0.10     2023-01-22 [1] CRAN (R 4.2.3)
##  RcppParallel        5.1.8      2024-07-06 [1] CRAN (R 4.2.3)
##  readr             * 2.1.4      2023-02-10 [1] CRAN (R 4.2.3)
##  readxl            * 1.4.3      2023-07-06 [1] CRAN (R 4.2.3)
##  rematch2            2.1.2      2020-05-01 [1] CRAN (R 4.2.3)
##  remotes             2.5.0      2024-03-17 [1] CRAN (R 4.2.3)
##  reshape2            1.4.4      2020-04-09 [1] CRAN (R 4.2.3)
##  respirometry      * 2.0.0      2024-07-18 [1] CRAN (R 4.2.3)
##  rlang               1.1.0      2023-03-14 [1] CRAN (R 4.2.3)
##  rmarkdown           2.21       2023-03-26 [1] CRAN (R 4.2.3)
##  robustbase          0.99-4-1   2024-09-27 [1] CRAN (R 4.2.3)
##  rstan             * 2.32.6     2024-03-05 [1] CRAN (R 4.2.3)
##  rstantools          2.4.0      2024-01-31 [1] CRAN (R 4.2.3)
##  rstudioapi          0.14       2022-08-22 [1] CRAN (R 4.2.3)
##  sandwich            3.1-0      2023-12-11 [1] CRAN (R 4.2.3)
##  sass                0.4.9      2024-03-15 [1] CRAN (R 4.2.3)
##  scales              1.3.0      2023-11-28 [1] CRAN (R 4.2.3)
##  seacarb             3.3.3      2024-02-15 [1] CRAN (R 4.2.3)
##  sessioninfo       * 1.2.2      2021-12-06 [1] CRAN (R 4.2.3)
##  shape               1.4.6.1    2024-02-23 [1] CRAN (R 4.2.3)
##  shiny             * 1.8.1.1    2024-04-02 [1] CRAN (R 4.2.3)
##  shinybusy         * 0.3.3      2024-03-09 [1] CRAN (R 4.2.3)
##  shinycssloaders   * 1.0.0      2020-07-28 [1] CRAN (R 4.2.3)
##  snakecase           0.11.1     2023-08-27 [1] CRAN (R 4.2.3)
##  SolveSAPHE          2.1.0      2021-05-01 [1] CRAN (R 4.2.0)
##  SRS               * 0.2.3      2022-03-27 [1] CRAN (R 4.2.3)
##  StanHeaders       * 2.32.9     2024-05-29 [1] CRAN (R 4.2.3)
##  stringi             1.7.12     2023-01-11 [1] CRAN (R 4.2.2)
##  stringr           * 1.5.1      2023-11-14 [1] CRAN (R 4.2.3)
##  survival            3.5-3      2023-02-12 [2] CRAN (R 4.2.3)
##  svUnit              1.0.6      2021-04-19 [1] CRAN (R 4.2.3)
##  tensorA             0.36.2.1   2023-12-13 [1] CRAN (R 4.2.3)
##  TH.data             1.1-2      2023-04-17 [1] CRAN (R 4.2.3)
##  tibble            * 3.2.1      2023-03-20 [1] CRAN (R 4.2.3)
##  tidybayes         * 3.0.6      2023-08-12 [1] CRAN (R 4.2.3)
##  tidygraph           1.3.0      2023-12-18 [1] CRAN (R 4.2.3)
##  tidyr             * 1.3.0      2023-01-24 [1] CRAN (R 4.2.3)
##  tidyselect          1.2.1      2024-03-11 [1] CRAN (R 4.2.3)
##  tidyverse         * 2.0.0      2023-02-22 [1] CRAN (R 4.2.3)
##  timechange          0.2.0      2023-01-11 [1] CRAN (R 4.2.3)
##  twosamples          2.0.1      2023-06-23 [1] CRAN (R 4.2.3)
##  tzdb                0.3.0      2022-03-28 [1] CRAN (R 4.2.3)
##  urlchecker          1.0.1      2021-11-30 [1] CRAN (R 4.2.3)
##  usethis           * 2.2.3      2024-02-19 [1] CRAN (R 4.2.3)
##  utf8                1.2.3      2023-01-31 [1] CRAN (R 4.2.3)
##  V8                  4.4.2      2024-02-15 [1] CRAN (R 4.2.3)
##  vctrs               0.6.1      2023-03-22 [1] CRAN (R 4.2.3)
##  vegan             * 2.6-6.1    2024-05-21 [1] CRAN (R 4.2.3)
##  viridisLite         0.4.2      2023-05-02 [1] CRAN (R 4.2.3)
##  withr               3.0.0      2024-01-16 [1] CRAN (R 4.2.3)
##  xfun                0.38       2023-03-24 [1] CRAN (R 4.2.3)
##  xml2                1.3.3      2021-11-30 [1] CRAN (R 4.2.3)
##  xtable              1.8-4      2019-04-21 [1] CRAN (R 4.2.3)
##  yaml                2.3.7      2023-01-23 [1] CRAN (R 4.2.3)
##  zoo                 1.8-12     2023-04-13 [1] CRAN (R 4.2.3)
## 
##  [1] C:/R
##  [2] C:/Program Files/R/R-4.2.3/library
## 
## ──────────────────────────────────────────────────────────────────────────────